Treatment planning

MatRad [23], an open-source treatment planning system (TPS), was developed for educational and research purposes. This also supports carbon-ion IMPT treatment planning. The energies of the carbon ions available in matRad are in the range of 115.23 MeV/u to 398.84 MeV/u corresponding to a penetration range of 32.68 mm to 294.25 mm in water. Furthermore, alpha-beta ratios of 10 and 2 were selected for the target volume and normal tissue, respectively. The dose grid was set to 3 mm × 3 mm × 3 mm in this study. The pencil beam scanning (PBS) method with spot scanning was used. Furthermore, PBS can realize conformal dose distributions for irregularly shaped tumors and improve the beam utilization [24].

Gaussian-shaped carbon-ion beam spot sizes range from 11.80 mm to 8.01 mm (full width at half maximum, FWHM) for energies in the range of 115.23 MeV/u to 398.84 MeV/u. A spot spacing of 3 mm in the lateral and longitudinal directions was employed. Carbon-ion IMPT treatment planning was performed in five patients with lung tumor lesions. For each patient, 4DCTs were acquired with 10 motion states. To reduce the consumption of computer memory, both 4DRO methods were developed using three significantly different breathing phases: 0% expiration phase, 80% expiration phase, and 40% inspiration phase. The 0% expiration phase was defined as the reference breathing phase.

For the PTV-based optimization, the PTV was generated by extending the union of the CTVs in all three phases with a 5- mm safety margin. The magnitude of the margins was derived from the original photon treatment planning [25]. The dose distributions were calculated for the three respiration phases without accounting for setup and range uncertainties. Both 4DRO methods were performed on the CTV for each phase without incorporating any margins. The setup uncertainties were assumed as ± 5 mm in ± x, ± y, and ± z directions. The range of the uncertainty parameters was set to ± 3.5% when the water-equivalent depth was computed. Each breathing phase had eight uncertain dose scenarios and one nominal dose scenario. The tumor size, location, magnitude of motion, and beam angle of the five patients under investigation are listed in Table 1. The CTV centroids in the right–left (RL), anterior–posterior (AP), and superior–inferior (SI) directions for the three breathing phases were recorded. The maximum centroid was used to determine the magnitude of the motion. The patient datasets included 4DCT images with CTV and OARs delineated by an experienced physician following the delineation guidelines provided in [26]. The prescription for all patients was 60 Gy(RBE)/15 fractions with threshold doses of 45, 20, and 30 Gy(RBE) for the spinal cord, lungs, and heart, respectively.The LEM was used to calculate RBE values.

4D robust optimization

Liu et al. [21] proposed a robust version of the 4D-optimization method and demonstrated that this approach can enhance the robustness of the cumulative 4D dose distribution. The physical absorbed dose distributions $\left(d_{i_k\mathrm{,}j}^s{w}_j^2\right)$ in different setups and range uncertainty scenarios were calculated for the breathing phases of 4DCT datasets, where wj represents the intensity weight of the beamlet j, and $d_{i_k\mathrm{,}j}^s$ represents the influence matrix describing the distribution of the beamlet j in the breathing phase k to voxel ik in the uncertainty scenarios. We define s as each uncertainty scenario considered during treatment planning optimization and define S as the set of all scenarios. During the optimization process, a non-negative vector $w_j^2$ was used to transform the constrained optimization problem (w≥0) into a non-constrained problem. The cumulative physically absorbed dose distributions were deformed to the reference breathing phase by using the registration algorithm in matRad. The cumulative RWD distributions $D_i^s$ were computed using the cumulative physically absorbed dose distributions. Specifically, $D_i^s$ were employed in a standard quadratic objective function to design a carbon-ion treatment plan. In the objective function, the worst-case scenario robust optimization mentioned above was conducted as follows: $\begin{array}{c}\text{min}f(w)={\displaystyle \sum _{N\in \text{Target}}\{q_{n\mathrm{,}\text{max}}\underset{s\in S}{{\displaystyle \max }}{\left(D_i^s-D_{\text{pre}}\right)}^2}\\ +q_{n\mathrm{,}\text{min}}\underset{s\in S}{{\displaystyle \min }}{\left(D_i^s-D_{\text{pre}}\right)}^2\}\\ +{\displaystyle \sum _{M\in \text{Oar}}{q}_{m\mathrm{,}\text{max}}}\underset{s\in S}{{\displaystyle \max }}H\left(D_i^s-D_{\text{o}}\right){\left(D_i^s-D_{\text{o}}\right)}^2,\end{array}$ (1) where D_pre and D_o denote the prescription dose in the target volume and dose constraints on the relevant tissues, respectively. Furthermore, N and M denote the numbers of target volumes and OARs, respectively. Parameters q_n,max, q_n,min, and q_m,max are the penalty factors for the target volumes and OARs. Furthermore, q_n,max=1.0, q_n,min=1.0, and q_m,max=1.0 were used. The step function H is mathematically defined as follows: H=1.0 when $D_i^s>D_{\text{o}}$ and H=0.0 otherwise. The ultimate objective of the optimization is to minimize the objective function f(w) to realize the optimal dose distributions over the relevant voxels.

The RP-4DRO method did not explicitly address the optimization of dose distributions in the individual breathing phases; it remained unclear whether robustness in individual respiratory phases can be ensured. Ge et al. [3] proposed the EP-4DRO method for solving the uncertainty problem in proton treatment planning. The RWD distributions $D_{i_k}^s$ in different setups and range uncertainty scenarios were calculated for all the breathing phases of 4DCT. The worst dose distributions for each breathing phase were calculated by selecting the worst objective score. Subsequently, the worst dose distributions for all the phases were simultaneously optimized. The objective function is expressed as follows: $\begin{array}{c}\text{min}f\left(w\right)={\displaystyle \sum _{N\in \text{Target}}({\displaystyle \sum _{k\in K}{q}_{n\mathrm{,}\text{max}}}}{\text{max}}_{s\in S}{\left(D_{i_k}^s-D_{\text{pre}}\right)}^2\\ {\displaystyle \sum _{k\in K}{q}_{n\mathrm{,}\text{min}}}{\text{min}}_{s\in S}{\left(D_{i_k}^s-D_{\text{pre}}\right)}^2\}+\\ {\displaystyle \sum _{M\in \text{Oar}}{\displaystyle \sum _{k\in K}{q}_{m\mathrm{,}\text{max}}}}{\text{max}}_{s\in S}H\left(D_{i_k}^s-D_{\text{o}}\right){\left(D_{i_k}^s-D_{\text{o}}\right)}^2.\end{array}$ (2) It should be noted that all the parameter settings were the same as those in the RP-4DRO method. During the optimization iteration, the objective function values for each dose scenario were calculated and the worst scenario was selected.

Plan evaluation

To maintain consistency, the resulting dose distributions for the volume of interests (VOIs) were evaluated and compared using the same parameters. Nine dose distributions for each breathing phase (6 setup uncertainty scenarios, 2 range uncertainty scenarios, and 1 nominal scenario) were calculated. The bands of dose-volume histograms (DVHs) derived from the 27 dose distributions (nine scenarios for three breathing phases) were analyzed to compare the robustness of the treatment plans. For PTV-based optimization, the dose distributions were recalculated for each breathing phase using the optimized set of spot intensities for each of the nine uncertainty scenarios. A narrower band in the DVH indicates greater robustness, implying that the dose distribution changes are relatively minimal under the influence of uncertainties. Evaluating DVHs with nine accumulated dose distributions is essential [3]. Thus, DVHs derived from the 9 accumulated dose distributions were also analyzed. Dose accumulation was assessed during the three breathing phases.

Dose distribution conformity index (CI) and target heterogeneity index (HI) were computed using the following equations: $CI=\frac{T{V}_{95\%}}{TV}\times \frac{T{V}_{95\%}}{V_{95\%}},$ (3) $HI=\frac{D_{5\%}-D_{95\%}}{D_{95\%}}.$ (4) Specifically, TV_95% denotes the target volume receiving a 95% prescription dose. Furthermore, V_95% represents the total volume receiving at least 95% of the prescribed dose; CTV is the target volume (TV) for the three methods. The CI is a measure of how well the target dose distributions conform to the target volume, with a value closer to 1.0 indicating greater conformity. HI reflects the homogeneity of the target dose distributions, with a value closer to 0.0 indicating greater homogeneity. To compare the target coverage, the index D_95% is used to represent the dose received by the 95% target volume. A value closer to the prescribed dose indicates better TV dose coverage of the target volume.

Normal tissue sparing was evaluated by comparing different treatment plans optimized using different methods. Evaluation indexes, such as lung D_mean, lung V20, lung V5, heart D_mean, heart V30, and spinal cord D_max were analyzed and compared among the aforementioned methods. Given that the dose distributions were non-Gaussian profiles, a statistical analysis was conducted using the Wilcoxon signed-rank test [27] to compare the PTV-based optimization and 4DRO methods. A statistically significant difference was considered when the P value was less than 0.05.

Exemplary patient

Figure 1 shows the DVH bands for the CTV and OARs derived from the 27 dose distributions for Patient 1 as an example. The results of the PTV-based optimization, RP-4DRO and EP-4DRO methods are shown in Fig. 1a, 1b, and 1c respectively. The dose distributions without considering the uncertainties for the reference phase (nominal scenario) are highlighted by black solid lines. For Patient 1, the bandwidths at D_95% for the PTV-based optimization, RP-4DRO, and EP-4DRO methods were 8.92, 8.77, and 5.29 Gy(RBE), respectively. Correspondingly, the bandwidths at D_5% were 2.11, 2.11, and 1.68 Gy(RBE). The planning objective of CTV coverage was set to $D_{95\%}/D_{\text{pre}}\ge 95\%$ in the nominal scenario. The CTV coverages in the nominal scenario for the PTV-based optimization, RP-4DRO, and EP-4DRO methods were 94.86%, 93.51%, and 96.73%, respectively.

Fig. 1Full-screen View

(Color online) DVHs of Patient 1 for the different treatment plans using the PTV-based optimization (a), RP-4DRO (b), and EP-4DRO (c) methods. The axial view of the accumulated 4D dose distribution overlaid on the reference phase is displayed on panels (d–f), corresponding to the different plans.

For the DVHs of the normal structures, both 4DRO methods realized narrower bands when compared to the PTV-based optimization in lung, heart, and spinal cord sparing. A nominal view of the accumulated 4D dose distributions in the reference phase for Patient 1 is presented in Fig. 1d, 1e, and 1f. The dose at the edges of the CTV for both 4DRO methods was reduced when compared to the PTV-based optimization. Notably, the RP-4DRO method (Fig. 1e) showed hot spots in the medial area of the CTV when compared with the EP-4DRO method (Fig. 1f).

Entire patient cohort

Figure 2　shows the target coverage indices for the nominal scenario.

Fig. 2Full-screen View

(Color online) Target coverage indices for the nominal scenario.

Compared with the PTV-based optimization, the target coverage indices realized by the RP-4DRO method were relatively low, except for Patient 3 (P=0.383). The target coverage indices were also below constraint ($D_{95\%}/D_{\text{pre}}\ge 95\%$). Conversely, the EP-4DRO method exhibited superior performance in terms of target coverage for all patients (P=0.225).

Figure 3 shows the robustness evaluation results and plan qualities with three optimization methods from five patients. The bandwidths at 95% and 5% CTV volumes for all three methods are shown in Fig. 3a and 3b, respectively. A lower value indicated a higher level of robustness. Based on the experimental findings, the dose distributions within the CTV derived from EP-4DRO exhibited narrower bandwidths than those achieved through PTV-based optimization, except for Patient 3. Similarly, the RP-4DRO method yielded narrower bandwidths in the dose distributions when compared to the PTV-based optimization, with the exception of Patients 2 and 3.

Fig. 3Full-screen View

(Color online) Bandwidth at 95% CTV volume (a) and bandwidth at 5% CTV volume (b) for the PTV-based optimization, RP-4DRO, and EP-4DRO approaches. The error-bar of heterogeneity (c) and conformity indices (d), respectively.

Figure 3c and 3d present the values of HI and CI for all 27 scenarios across the five patients. Error bars have been added to the figure to indicate one standard deviation of the results. For the EP-4DRO method, the HI values were generally lower than those for the RP-4DRO method and PTV-based optimization (P=0.345), except for Patient 4. Lower values indicate better dose homogeneity at the target boundaries. Conversely, the RP-4DRO method yielded relatively higher HI values, indicating worse dose homogeneity in the individual respiration phases (P=0.107). As shown in Fig. 3d, the CI values obtained using the EP-4DRO method are higher compared to the PTV-based optimization (P=0.160) and RP-4DRO method. The CI values produced by the EP-4DRO method were higher than those of the PTV-based optimization except for Patient 4 (P=0.005).

As shown in Fig. 4, the dose-volume indices for the OARs are also calculated. A comparative analysis revealed that the RP-4DRO method delivered lower doses to the lung ((left/right) lung D_mean, P<0.05; left lung V20, P=0.068; (left/right) lung V5, P<0.05; and right lung V20, P<0.05). As shown in Fig. 5, both 4DRO methods yielded lower indices for the heart and spinal cord when compared to the PTV-based optimization (P<0.05), indicating the heart and spinal cord received less doses. Notably, for the five patients, the RP-4DRO method consistently exhibited the lowest doses to the lungs and heart across the individual breathing phases. However, the RP-4DRO method did not demonstrate superiority in maintaining lower doses delivered to the spinal cord when compared with the EP-4DRO method. The P values indicate that the differences in the OARs indices between the two 4DRO methods were statistically significant.

Fig. 4Full-screen View

(Color online) Dose distributions in lung for the PTV-based optimization, RP-4DRO, and EP-4DRO methods.

Fig. 5Full-screen View

(Color online) Dose distributions in heart and spinal cord for the PTV-based optimization, RP-4DRO, and EP-4DRO methods.

The aforementioned results illustrate 27 dose distributions for the three respiratory phases. However, for further analysis in a clinical implementation, it is necessary to compare and calculate the accumulated dose distributions over the three respiratory phases.

The bandwidth at 95% CTV volume and target coverage of the cumulative distributions are shown in Fig. 6. The RP-4DRO method showed a smaller bandwidth at 95% CTV volume when compared to the PTV-based optimization (P=0.022). Conversely, the EP-4DRO method does not guarantee small bandwidths for Patient 2 and Patient 5 when compared with the PTV-based optimization (P=0.272). The target coverages generated by the RP-4DRO method (P=0.45) and EP-4DRO method (P=0.066) were lower than that of the PTV-based optimization. The dose distributions in OARs are illustrated in Fig. 7. Compared with the PTV-based optimization, both 4DRO methods achieved lower doses in the lung, heart, and spinal cord (P<0.05). Notably, the EP-4DRO method resulted in the lowest doses in the lungs and heart. In contrast, the EP-4DRO method was not superior in maintaining the lowest cumulative dose distribution delivered to the spinal cord. This implied that the accumulated dose distributions generated by both 4DRO methods are more conservative.

Fig. 6Full-screen View

(Color online) Accumulated results of robustness and target coverage for the PTV-based optimization, RP-4DRO, and EP-4DRO methods.

Fig. 7Full-screen View

(Color online) Accumulated dose distributions in lung, heart, and spinal cord for the PTV-based optimization, RP-4DRO, and EP-4DRO methods.