Beam‐specific planning target volumes incorporating 4D CT for pencil beam scanning proton therapy of thoracic tumors

The purpose of this study is to determine whether organ sparing and target coverage can be simultaneously maintained for pencil beam scanning (PBS) proton therapy treatment of thoracic tumors in the presence of motion, stopping power uncertainties, and patient setup variations. Ten consecutive patients that were previously treated with proton therapy to 66.6/1.8 Gy (RBE) using double scattering (DS) were replanned with PBS. Minimum and maximum intensity images from 4D CT were used to introduce flexible smearing in the determination of the beam specific PTV (BSPTV). Datasets from eight 4D CT phases, using ±3% uncertainty in stopping power and ±3 mm uncertainty in patient setup in each direction, were used to create 8×12×10=960 PBS plans for the evaluation of 10 patients. Plans were normalized to provide identical coverage between DS and PBS. The average lung V20, V5, and mean doses were reduced from 29.0%, 35.0%, and 16.4 Gy with DS to 24.6%, 30.6%, and 14.1 Gy with PBS, respectively. The average heart V30 and V45 were reduced from 10.4% and 7.5% in DS to 8.1% and 5.4% for PBS, respectively. Furthermore, the maximum spinal cord, esophagus, and heart doses were decreased from 37.1 Gy, 71.7 Gy, and 69.2 Gy with DS to 31.3 Gy, 67.9 Gy, and 64.6 Gy with PBS. The conformity index (CI), homogeneity index (HI), and global maximal dose were improved from 3.2, 0.08, 77.4 Gy with DS to 2.8, 0.04, and 72.1 Gy with PBS. All differences are statistically significant, with p‐values <0.05, with the exception of the heart V45 (p=0.146). PBS with BSPTV achieves better organ sparing and improves target coverage using a repainting method for the treatment of thoracic tumors. Incorporating motion‐related uncertainties is essential. PACS number: 87.55.D

uses a perpendicular search radius, typically fixed at 5 mm, over the entire beam path to add treatment margin to the target along the beam path and was implemented in the treatment of thoracic tumors by Moyers et al. (2) using the DS technique. In contrast to the implicit smearing embedded in the compensator design and direct manipulation of the distal and proximal ranges of DS beams, Park et al. (3) introduced a beam-specific PTV (BSPTV), initially proposed by Rietzel and Bert, (4) to explicitly include variation of water-equivalent path length (WEPL) along each beam direction; in this manner the BSPTV can be used in pencil beam scanning (PBS) planning single field optimization (SFO). In the BSPTV method, (3,4) distal and proximal water-equivalent treatment margins (WETM) are converted to geometric treatment margins (GTM) that are calculated according to local tissue heterogeneity and added beyond the target to achieve a smearing effect in PBS, thus accounting for WET variations related to the fixed value of misalignment of tissue from motion and setup.
Neither Moyer's nor Park's approaches account for the fact that the magnitude of motion varies within different anatomic regions. The modeling of treatment margin due to organ motion can be further improved by using maximum intensity projections (MIP) and minimum intensity projections (MinIP) obtained from 4D CT. (5) Matney et al. (6) and Flampouri et al. (7) used 4D CT for evaluation of robustness and for construction of a BSPTV in DS treatment plans, respectively.
Respiratory motion patterns change among different segments of the thorax region. (8) The magnitude of motion at a patient's surface (i.e., at the beam entrance) can be significantly smaller than that of tumor motion. (9,10) Therefore, a fixed smearing value based on the magnitude of target motion would typically be too large for proximal regions along the beam path, and a BSPTV using a fixed smearing radius corresponding to tumor motion magnitude would be larger than needed. Multiple researchers have attempted to refine the BSPTV method using 4D CT for PBS treatment. Graeff et al. (11) attempted to extend BSPTV to multiple field optimization (MFO) IMPT fields in GSI's in-house TPS. Knopf et al. (12) calculated BSPTV as the union of multiple treatment targets (CTV + margins) over different phases of 4D CT using the PSI and NIRS in-house deformable registration algorithms.
PBS treatment planning for thoracic tumors can achieve better sparing of organs at risk (OAR) than IMRT techniques. (13,14,15) Despite the dosimetric advantages of PBS, DS and IMRT remain the methods of choice for treating complex thoracic tumors. This is primarily because commercial treatment planning systems (TPS) are unable to: 1) calculate proton spot delivery interplay with organ motion, which can lead to overdose and underdose within the treatment volume, and 2) determine appropriate treatment margins around the treatment target due to uncertainties of motion, stopping power, and setup, which might result in partial anatomical miss or underdose of the target volume. Although Li et al. (16) recently evaluated the systemic errors in 3D dose calculation that appear in the PBS treatment of thoracic tumors, there is no simultaneous comparisons of OAR sparing and plan robustness between PBS and DS treatments in this disease site. (17) In this work, we describe the development of a patient field-specific BSPTV that incorporates respiratory motion and stopping power uncertainties using 4D CT for PBS treatment of thoracic tumors. As reported by previous researchers, (11,12) we calculated the change of waterequivalence path length (WEPL) associated with organ motion and converted such WETM to GTM based on local tissue densities. In contrast to the conceptual demonstration by Graeff et al. (11) for MFO application, we focus on SFO application in multiple thoracic patients. And unlike Knopf et al., (12) who conceptualize the method of using the union of (multiple deformable targets + treatment margins) over different breathing phases in two different in-house TPS of PSI and NIRS, our treatment margin is applied to average CT images and the iCTV, which is the union of multiple CTV targets rigidly registered from eight phases in a commercially available TPS (Varian Eclipse, Palo Alto, CA). Compared to Knopf's and Graeff's approaches, we need only to calculate WETM and convert WETM to GTM once instead of eight or ten times (the # of breathing phases). This may lead to slightly larger target volumes than those of Knopf and Graeff, similar to the more generous margin calculation method proposed by Flampouri et al. (7) for DS treatment. To illustrate the potential advantages of BSPTV using 4D CT for the treatment of thoracic tumors, we will demonstrate how such a method is equivalent to "flexible smearing" in DS and can generate a more appropriate treatment target volume for a typical patient.
A number of pioneering researchers have demonstrated the magnitude of overdose and underdose within the target volume caused by PBS delivery interplay with organ motion and how such deviation from the prescribed dose distribution is washed out with multiple fractions. (18,19,20,21) In this study we demonstrate such an interplay phenomenon for a representative thoracic patient.
The focus of this study is to establish a method of BSPTV calculation and organ interplay so that we can systemically report the dosimetric difference of PBS and DS treatment methods for 10 thoracic tumor patients. Such methods and dosimetric studies can also serve as the platform and baseline for future improvements that can be brought by various motion mitigation strategies, such as multiple beam paintings, and gating and tracking techniques. (22)

II. MATERIALS AND METHODS
In conventional radiation therapy, the internal clinical target volume (iCTV) is delineated by a radiation oncologist and by the union of CTVs on the corresponding 8-or 10-phase images of a 4D CT on average CT images. To account for random and systematic errors of patient setup, an additional margin of 5-8 mm is used to expand an iCTV to a PTV.
Our methodology to calculate WETM and convert WETM to GTM is based on the phase and average images of 4D CT and implemented as a standalone MATLAB (MathWorks, Natick, MA) program. The generated BSPTV were imported into Eclipse TPS (Varian Medical Systems) for subsequent optimization, thus providing better PBS treatment plan of thoracic tumors. Instead of using a "fixed smearing" distance to search for maximal WEPL change over a beam path, the maximal difference of WEPL between 8 4D CT phases and the average CT images was used as WETM along each beam path. The distal and proximal water-equivalent margins (DM we and PM we ) due to organ motion (WETM) were calculated as the WEPL difference between the MIP, MinIP, and the average of the eight phases of the 4D CT, respectively (Eq. (1) and (2)). In this manner, different smearing radii along the beam path (i.e., "flexible smearing") were effectively used to calculate WETM. Such WETM were later implemented as additional voxels beyond treatment targets (i.e., GTM) based on the local Hounsfield units in the average CT images and the projected voxel size over the beam path, whose summation should equal to WEPL differences (Eq. (3) and (4)). The steps of conversion of WETM to GTM were omitted in the original smearing approaches by Urie et al. (1) and Moyers et al., (2) but implemented in recent works of others. (3,4,7,11) Figure 1(a) shows an example of how a fixed smearing of 5 mm radius over average CT images overestimates a BSPTV compared to that of flexible smearing over 4D CT images. Two BSPTVs, generated using fixed smearing and flexible smearing, are shown for a gantry angle of 270°. Because the motion of the ribs is less than that of the target, a fixed smearing of 5 mm radius overestimates the distal and proximal BSPTV margins. For the regions close to the heart, however, the distal BSPTV margin from flexible smearing can be slightly larger than that from fixed smearing because cardiac motion can exceeds 5 mm. In addition, cardiac motion is not inadequately modeled by breathing phase-based 4D CT images. Figure 1(b) demonstrates how the beam-specific range uncertainties are broken down into three aspects: internal motion (Motion), stopping power ratio (SPR) of medium-to-water versus Hounsfield units in CT images, and the accuracy of image guidance (Setup). The BSPTV calculation uses 3% of the proton range plus 1 mm related to CT/stopping powers and, for patient setup variation, 2 mm along AP/RL directions and 3 mm along the SI direction. As these three sources of uncertainties are most likely independent, a quadratic method was primarily used for the summation of three uncertainties; however, a linear summation was also used to ensure  target coverage, providing there was no overlap with adjacent organs at risk. Figure 1(c) shows the contributions to the BSPTV as a function of gantry angle. The contribution from Motion is larger than from CT/SPR or Setup, and contributes over 90% of the Sum as determined by the following equations: where i are the voxels along the beamlet from the patient surface to the proximal and distal surfaces of iCTV in Eqs. (1) and (2), respectively; j are the amount of voxels to be determined in the distal and proximal GTM regions in Eqs. (3) and (4), respectively; and where l and dl are voxels and projected voxel dimension along the beam path in average CT, respectively.
We further define two overlap volumes -one between the BSPTV and organs, and the other between the beam path outside of the BSPTV and organs. While the beam path overlap volume may receive less than the full prescription dose of 66.6 Gy, minimizing the overlapped volumes along with the BSPTV will allow for sparing of lung, heart, and cord prior to the step of plan optimization. Figure 1(d) shows how the MATLAB code generates an optimized BSPTV and overlapped volume at an ideal gantry angle. In this figure, a gantry angle of 300° was selected by both MATLAB and the dosimetrist, since it corresponded to the long axis of the target and optimal BSPTV and overlap volumes. In contrast, a gantry angle of 180° chosen by the dosimetrist is suboptimal, since the BSPTV and overlap volumes are both relatively large.
Motion, CT/SPR, Setup, and Quadratic and Linear Summations BSPTVs were used in the optimization step of PBS treatment planning for each of the 10 patients evaluated. The plans were optimized with goals of covering 99% of the BSPTV with 99% of the prescription dose of 66.6 Gy in 1.8 Gy (RBE) per fraction, though the coverage of BSPTV could be selectively compromised if there was a critical OAR constraint. The OAR constraints within the optimizer were set according to the constraints used in the RTOG 1308 protocol, (23) a phase III randomized trial comparing photon versus proton chemoradiotherapy for inoperable stage II-IIIB NSCLC. Two fields were used in each PBS treatment plan, and each was individually optimized and summed at the last step without any further optimization. The treatment angles were kept the same in PBS plans as they were in DS plans unless the plan optimizer preferred other angles.
Any comparison of treatment plans depends on the characteristics of beam lines used. The spot and Bragg peak characteristics of the IBA pencil beam scanning beam line have been reported by Farr et al., (24) Grevellot et al., (25) and Lin et al., (26)(27) while the source size and Bragg peaks of the double scattering technique have been described by Slopsema et al. (28) The spot sigma at most patient surfaces decrease from 7 to 3 mm for proton energies from 100 to 225 MeV. In order to treat tumors shallower than 75 mm (i.e., the proton range of 100 MeV beams), we utilize a U-shaped bolus (Both et al. (29) ) with a WET of 75 mm placed on the top and lateral sides of patient. Our proton couch also has a built-in bolus with a WET thickness of 65 (or 75 mm with an overlay) to enable treatment of superficial posterior tumors. The dose algorithm PCS13.0.24 was used for all PBS and DS calculations of Eclipse TPS. The commissioning of Eclipse has previously been reported by Zhu et al. (30) To facilitate plan comparison between PBS and DS techniques, several dosimetric indices of target iCTV were used in this work, as described by Flampouri et al.: (7) the volume of 100% isodose line, IDL 100% ; the Heterogeneity index, HI = (D 5% -D 95% )/D prescription ; the Conformity index, CI = IDL 100% /V ctv ; Coverage quality, CQ = D min /D prescription ; and Integral dose, ID = V body × D avg .
Lomax et al. (31) displayed uncertainty bands on the DVH curve to designate the uncertainties due to stopping power and patient setup. The uncertainty band method was also implemented by Lin et al. (32) and Vargus et al. (33) for OAR and target coverage for patient setup variation. In this work, a total of 8 × 12 × 10 = 960 plans were evaluated for 10 patients, comprising permutations of 8 4D CT phases, ± 3% stopping power, and ± 3 mm x/y/z in patient setup. Percentiles of 25%-75% of all the permutations from the average DVH of OARs and iCTV were calculated for a representative patient and all 10 patients (Fig. 2 and Fig. 3). The difference between PBS and DS was also calculated and displayed.
Ten consecutive patients treated with DS for stage III locally advanced non-small cell lung cancer were selected for the evaluation of PBS plans. All patients had mediastinal nodal metastases with various primary tumor locations. The comparisons of DVH criteria were evaluated using a paired t-test for statistical significance.
The method we used to evaluate beam interplay is similar to Zou et al. (33) for conventional lung SBRT applications and other researchers for PBS applications. (19,18,11,12,21) A representative patient is shown in Fig. 4, with a breathing period of ~ 3.5 s and delivery durations of ~ 60 s and ~ 46 s for the two beams. Switching time between energy layers, slew duration between  spots and spot duration of the same energy layer were extracted from the beam delivery log files. For different treatment fractions or beam paintings within the same fraction, the beams start from a random position in the 3.5 s breathing period. PBS spots were, therefore, grouped into the 8 different breathing phases. A single treatment plan was, therefore, split into eight plans outside of Eclipse using our in-house MATLAB programs and a DICOM editor provided by Ion Beam Applications (Louvain-La-Neuve, Belgium). These plans were subsequently reimported into Eclipse TPS to calculate the dose distribution of each CT phase. The eight dose distributions were deformed and registered to the exhale-50 phase and summed in MIM Maestro (MIM Software Inc., Cleveland, OH). Table 1 shows the magnitude of motion at three locations along the beam path for the 10 consecutive patients evaluated. Tumor motion was larger than that of the ribs or at the patients' surface (beam entrance), and it was largest in the SI direction (~ 8 mm) followed by AP (~ 4 mm) and RL (~ 2 mm) directions. In contrast, the average motion of the ribs and patient surface at the beam entrance were less than 1.5 mm and 1 mm, respectively. Out of 10 patients, all of whom were treated with two-field plans, six plans were treated using posterior and posterior oblique beams only, two plans involve a lateral beam, and the other two plans included an anterior beam. With supine positioning, there is minimal motion of ribs and posterior beam entrance. In contrast, when lateral or anterior beams were used, the motion was larger. Table 2 shows BSPTV as a function of source of uncertainty. In the summed BSPTV, motion-induced range uncertainties were the largest contribution to treatment margins. Figure 2 shows a representative patient treated with DS who would have benefited from a PBS technique. PBS spared both the high-and low-dose regions and allows reduced doses to the lung, heart, cord, and esophagus compared with DS technique, while still maintaining 99% dose coverage of 99% of the iCTV. Furthermore, dose inhomogeneity and hotspots within the target were minimized, whereas target coverage was maintained for both PBS and DS (within ~ 2% dose) when 4D CT, stopping power, and patient setup uncertainties were incorporated. Figure 3(a) shows the average, 25th, and 75th percentile DVH bars for lung, heart, cord, and esophagus of all 10 patients for PBS (red) and DS (blue). P-values > 0.05 suggest that the PBS advantage observed may not be statistically representative of each individual patient. PBS plans were superior to DS in all cases. The greatest improvements with PBS was primarily within the low (V5) to moderate (V10 and V20) dose regions, with a reduction in mean lung dose (MLD) by 2.3 Gy from 16.4 Gy with DS to 14.1 Gy, and heart V30 from 10.4% to 8.1% Table 1. Magnitude of motion in mm observed along the beam path for the target, rib, and entrance. with PBS (Table 2). There is a large variation in heart dose among the 10 patients. Because one patient had almost no dose to the heart, the maximal heart dose average over 10 patients fell outside of the 25 and 75 percentiles for both PBS and DS treatments. Excluding this patient, PBS had a maximal heart dose of 70.9 Gy averaged over nine patients, with 25 and 75 percentiles from 70.3 to 71.3 Gy, while DS had a maximal heart dose of 74.1 Gy averaged over nine patients, with 25 and 75 percentiles from 70.8 to 75.1 Gy. The reduction in heart V45, was not statistically significant (p = 0.146).

III. RESULTS
In keeping with the dose constraints from RTOG 1308, only the maximum dose was constrained for esophagus. Despite this, PBS still resulted in a reduction in the average volume of esophagus irradiated. A p-value of 0.049 indicates that the D max for PBS can be larger than for DS for some patients, although this occurs only when D max is well below 74 Gy. Because one patient had a maximal esophagus dose of 48 Gy in PBS treatment, the average over 10 patients fell outside of the 25 and 75 percentiles. Excluding this patient, PBS had a maximal esophagus dose of 70.0 Gy averaged over nine patients, with 25 and 75 percentiles from 69.5 to 71.2 Gy, while DS had a maximal esophagus dose of 72.6 Gy, with 25 and 75 percentiles from 70.5 to 73.5 Gy. The reduction in maximum heart and esophageal dose for certain patients from 72 Gy to 68 Gy (Table 2 and Fig. 3) suggests that additional patients would benefit from PBS. Table 3 shows that the integral dose (ID) is reduced by nearly 17% in PBS compared with DS. Additionally, the irradiated volume of the prescription dose (CI) is only reduced by 10%.
The target coverage quality CQ is identical between PBS and DS, as 99% of the iCTV was required to receive the prescription dose with both techniques. In contrast, PBS results in a significantly more homogenous dose, with a homogeneity index HI of 0.04 compared with 0.08 for DS. Except for the marginal advantage of CI (p = 0.056), all the differences above were highly statistically significant, with p < 0.01.
The robustness of iCTV coverage is shown in Fig. 3(b) for all 10 patients. For all PBS and DS plans, the iCTV receives no less than 97% of the prescription dose when all permutations of motion, CT/stopping power, and patient setup uncertainties are included. Despite many clinical fears over overdose/underdose within the treatment target and target miss, it seems that PBS can potentially provide robust target coverage comparable or superior to that of DS when uncertainties are properly accounted for in the BSPTV with a repainting method discussed below.

IV. DISCUSSION
Flampouri et al. (7) proposed two methods to improve BSPTV in DS over the approach of Park et al. (3) by incorporating respiratory motion. The first method uses MIP and MinIP data from 4D CT to determine the distal and proximal margins, respectively. In the second, the 4D CT phase that generates the largest margins is selected for each beamlet (beam path) to the target. The beamlet method results in a smaller BSPTV than the MIP/MinIP method. Unlike DS, however, in which all energy layers are delivered within 0.1 s, PBS delivers different energy layers sequentially over minutes. Thus delivery of the distal and proximal layers are out of phase and, in contrast with DS, more generous margins must be determined from MIP/MinIP instead of the worst 4D CT phase in PBS. We demonstrated that the contribution from motion is the dominant effect over stopping power and setup uncertainties to the treatment margins in PBS and concluded that "flexible smearing" based on MIP/MinIP is crucial to the success of PBS plan optimization. We further included the overlapped volumes of beam path with organs. Unlike DS, which always over treats proximal organs, sparing of such OARs by PBS beam can be substantial. Our finding that the iCTV is more homogenously covered in PBS than in DS does not account for organ interplay with beam delivery in PBS. Accounting for interplay, the target dose distribution in each of 37 fractions of PBS is more heterogeneous than the planned distribution. The breathing period of ~ 3.5 s is comparable with the switching time between energy layers (1-2 s for our PBS system) and also comparable with the time required to paint each energy layer (with 3-5 ms per spot and 2-4 ms slew duration between adjacent spots with a total of ~ 6000 spots over ~ 20 energy layers). Our study indicates that an energy layer can be delivered in 1-4 out of 8 different breathing CT phases of different breathing cycles. As our spot size is relatively small, more fractions might be needed to achieve total dose homogeneity compared with centers employing larger spot sizes. Figure 4 shows the interplay on an axial slice; underdose (95%) and overdose (110%) are present for a single fraction without any motion mitigation. As a worst case scenario, dose heterogeneities from 90%-120% can occur in a single fraction. However, such a magnitude of overdose and underdose within the target will be reduced when more treatment fractions are used. With 4 fractions, the extent of such overdose drops to 105%, with no 95% underdose, and with 8 fractions the 105% overdose volume decreases even further. The summated distribution from the deformable registered dose distributions of 4 fractions accounting interplay has a slightly larger area of 105% isodose and slightly smaller minimal dose than the planned dose distribution. With 8 fractions, there was no observable difference for the DVH of the iCTV that accounted for interplay and the original treatment plan without interplay. As multiple beam paintings of four are clinically achievable, we believe that the dosimetric advantage of PBS treatment in this patient will be able to be maintained despite the interplay phenomenon. These observations are consistent with that reported by many researchers. (4,19,18,21,20) Because the breathing period and beam duration can vary from patient to patient, detailed assessment of the interplay effect for multiple patients and multiple motion mitigation strategies is the subject of current investigation within our group. Interplay is not an issue in DS, since all energy layers are painted every 0.1 s.
Any plan comparison between PBS and DS techniques depends on the spot size of the PBS system, the placement of bolus/range shifter, the source size of the DS system, and the choice of prescription dose. The in-air proton spot size of our beam lines is smaller than those of the Hitachi system reported by Gillin et al. (35) Furthermore, because the built-in posterior bolus and U-shaped anterior and lateral bolus are closer to the patient surface than range shifters mounted on the PBS nozzle, the spot size incident on the patient is smaller than systems that use nozzle-mounted range shifters. Therefore, conclusions drawn from this study may not be applicable to other beam lines.

V. CONCLUSIONS
The use of a beam-specific PTV based on respiratory motion along the beam path derived from 4D CT and incorporating CT/stopping power and patient setup uncertainties can achieve better organ sparing with PBS compared to DS. Additionally, our study indicates that PBS planning based on BSPTV method can potentially achieve both better target homogeneity and robustness of target coverage using a repainting delivery method. Finally, beam angle optimization provides a reliable method to minimize doses to the critical organs.