Four‐dimensional computed tomography‐based ventilation imaging in intensity‐modulated radiation therapy treatment planning for pulmonary functional avoidance

Abstract Purpose To incorporate four‐dimensional computed tomography (4DCT)‐based ventilation imaging into intensity‐modulated radiation therapy (IMRT) treatment planning for pulmonary functional avoidance. Methods and Materials Nineteen locally advanced lung cancer patients are retrospectively studied. 4DCT images are employed to create ventilation maps for each patient via a density‐change‐based algorithm with mass correction. The regional ventilation is directly incorporated into the mathematical formulation of a direct aperture optimization model in IMRT treatment planning to achieve functional avoidance and a voxel‐based treatment plan. The proposed functional avoidance planning and voxel‐based planning are compared to the conventional treatment planning approach purely based on the anatomy of patients. Paired sample t‐tests are conducted to see whether dosimetric differences among the three approaches are significant. Results Similar planning target volume (PTV) coverage is achieved by anatomical, functional avoidance, and voxel‐based approaches. The voxel‐based treatment planning performs better than both functional avoidance and anatomical planning to the lung. For a total lung, the average volume reductions in a functional avoidance plan from an anatomical plan, a voxel‐based plan from an anatomical plan, and a voxel‐based plan from a functional avoidance plan are 7.0%, 16.8%, and 10.6%, respectively for V40; and 0.4%, 6.4%, and 6.0%, respectively for mean Lung Dose (MLD). For a functional lung, the reductions are 8.8%, 17.2%, and 9.2%, respectively, for fV40; and 1.1%, 6.2%, and 5.2%, respectively, for functional mean lung dose (fMLD). These reductions are obtained without significantly increasing doses to other organs‐at‐risk. All the pairwise treatment planning comparisons for both total lung and functional lung are statistically significant (p‐value α=0.05. From these results, we can conclude that voxel‐based treatment planning outperforms both anatomical and functional‐avoidance planning. Conclusions We propose a treatment planning framework that directly utilizes functional images and compares voxel‐based treatment planning with functional avoidance and anatomical treatment planning.


INTRODUCTION
Pulmonary ventilation imaging provides useful information on regional lung function. An emerging technique to compute lung ventilation map has been proposed [1,2] that uses four-dimensional computed tomography (4DCT) image datasets. Because 4DCT images are routinely acquired for lung cancer patients in current clinical practice, ventilation computation involves mere imaging processing. Compared to other imaging modalities, such as single photon emission computed tomography (SPECT) [3] and magnetic resonance (MR) [4] , 4DCT ventilation imaging has a higher resolution, lower cost, shorter scan time, and greater availability. [5] Research on 4DCT-based ventilation imaging has focused on its methodology, [6] validation and clinical usage. Several recent validation studies [5,7,8] have compared 4DCT-based ventilation with both regional nuclear medicine ventilation-perfusion imaging and the global pulmonary function test, and generally showed good agreement with the different techniques in assessing lung function. With respect to clinical usage, Vinogradskiy et al. [9] uses 4DCT-based ventilation maps to evaluate lung function changes throughout the course of radiotherapy treatment, and Vinogradskiy et al. [10] incorporates 4DCT ventilation-based function to better correlate lung dose with radiation pneumonitis. Another clinical use of 4DCT-based ventilation map was to create functional avoidance plans. [11,12] These studies generally involved generating a new contour corresponding to a highly functional lung region based on the ventilation map and imposing new dose-volume constraints on the contour so that desired doses would be delivered to the region. Although this approach did take into account heterogeneous lung function, it took a binary view of the entire function spectrum, essentially assuming function uniformity in the high-functioning lung as well as the low-functioning lung.
Our proposed methodology in this study instead takes full advantage of the entire function spectrum provided by the ventilation map and directly incorporates lung ventilation at every voxel into a mathematical model for inverse planning. We use regional ventilation to modulate the penalty scalars associated with lung overdose on a voxel-by-voxel basis so that the voxels with high function are given relatively larger penalties than those with low function. This could thus potentially steer dose away from well-ventilated lung regions.
We remark that the proposed methodology is similar to that of St-Hilaire et al. [13] But instead of using SPECT perfusion information in an aperture-based IMRT model, we consider 4DCT-based ventilation in a beamlet-based model.

Patient CT and contours
Nineteen locally advanced lung cancer patients treated at University of Maryland School of Medicine are retrospectively studied. Each of these patients has a free-breathing CT and 4DCT scan prior to treatment. The gross tumor volume (GTV), clinical target volume (CTV), internal target volume (ITV), planning target volume (PTV), and organs-at-risk (OARs), including heart, esophagus, spinal cord, and total lung (minus GTV) are contoured on the planning CT image. The GTV to CTV margin is 6 mm. The ITV is obtained by adding a margin of 6 mm isotropically to the CTV, and the PTV is formed by expanding the ITV isotropically by a margin of 5 mm.

4DCT-based ventilation imaging
Using the 4DCT image data set for each patient, a ventilation map is calculated via a density-change-based algorithm with mass correction. [2] Lungs are segmented on the CT images at the inhale and exhale phases of the respiratory cycle. The images at these two phases are matched to each other using a deformable image registration technique. Let V IN and V EX be the inhale and exhale volumes for the voxel, and IN and EX are the voxel densities at the inhale and exhale phases. Let parameter k be the ratio of the entire lung mass between the inhale and exhale phases. It is introduced in Equation (1) to account for the tissue-based mass changes that occur over the respiratory cycle.The local ventilation for a given voxel is estimated by A lung can be treated as a composite of mainly air and tissue. Let V VOX , V AIR , and V TISSUE be the volumes for the voxel and its air and tissue components, respectively. The density of a voxel is approximated by where AIR ≈ 0.001g∕cm 3 and TISSUE ≈ 1.1g∕cm 3 .

Monte-Carlo dose calculation
An in-house Monte Carlo Kernel-based Superposition (MCKS) dose calculation engine is applied to compute beamlet dose distributions for each beam orientation. The dose deposited to every voxel in each structure (PTV or OAR) by each beamlet (namely the dose deposition matrix) is obtained as input to the mathematical model for inverse planning. The voxel size is 4 mm × 4 mm × 2.5 mm. The same beam orientations that are used clinically to treat the patients are used in this study, which ranged from 4 to 10. The beamlet size is 1 cm by 1 cm. Depending on the anatomy of the patient, a field size of 12 cm × 12 cm, 16 cm × 16 cm or 20 cm × 20 cm is calculated to ensure adequate coverage of the tumor for each patient. Thus the total number of beamlets is between 576 and 2560.

Mathematical model
Let S be the set of all structures in the treatment volume. Each treatment volume consists of target structures (PTV),organs-at-risk (OARs),critical structures,and normal tissues. Esophagus, heart, and lung are examples of OARs, while the brain and spinal cord are in the critical structures category. Let T ⊆ S be the set of target structures, O ⊆ S be the set of OARs, and C ⊆ S be the set of critical structures. Let V s be the set of voxels for each s ∈ S. Let B be the set of beam angles. Each beam consists of m rows and n columns and forms a rectangular grid of beamlets. A two-dimensional matrix, which is called intensity matrix, is often used to specify the position and intensity of all beamlets. Each element of this intensity matrix is called a bixel. Let N be the set of all bixels where N = B × {1, .., m} × {1, .., n}. let K l be the set of allowable apertures for each beam l ∈ B, let K ≡ ∪ l∈B K l be the set of all apertures, and let A k ⊂ N be the set of beamlets that are exposed in aperture k ∈ K. Let parameter D ils be the dose deposition matrix for all i ∈ A k , l ∈ V s , s ∈ S. Let decision variables a k and z ls be the weight of aperture k ∈ K, and dose received by each voxel l ∈ V s and s ∈ S, respectively. z ls is defined as follows : Similar to this research, the optimization model in Romeijn et al. [14] describes a column generation approach to radiation therapy treatment planning using aperture modulation. They used head-and-neck cancer patients data to study the quality of the treatment plans generated using this optimization approach. However, in this research we use lung cancer patients data. The weighted optimization model (hereafter restricted problem) is shown below: Objective (4) is the weighted sum of target structures, OARs, and critical organs. For target structures, our goal is to minimize the maximum variation from the prescribed dose pd s . We want to minimize the mean dose for OARs. For critical organs, we wish to minimize the max dose.The rest of the volume of the treatment region is the normal tissue. We wish to minimize the maximum dose for normal tissue. Constraint set (5) guarantees that s percentage of target structure receives the prescribed dose pd s , and s is a free variable for all s ∈ T. [14] Constraint set (6) ensures the non-negativity of aperture weights.
The number of apertures produced in the restricted problem is much too large. However, it is expected that only few apertures will actually be used. [14] In this research, we start solving restricted problem by choosing a limited set of apertures. Then to check the optimality of the restricted problem, a subproblem called the pricing problem is solved to try to identify the apertures to enter into the restricted problem. [15] Let binary variable x i represents an aperture such that 1, if beamlet i ∈ N is exposed in the aperture, 0, otherwise.
Let X be the set of all binary vectors x that represents an allowable aperture. Then the pricing problem is defined as follows: If the value of the above pricing problem (7) is nonnegative, then the current solution is optimal. Otherwise, we include the generated new aperture to the restricted problem. We then re-optimize the new restricted problem and repeat it until there are no such apertures left that can improve the solution of the restricted problem.

Plan generation
The clinical dose-volume goals for the target volume and OARs are listed in Table 1. These goals are applied in anatomical (hereafter Anatom), functional avoidance (hereafter FuncAvoid), and voxel-based (hereafter Voxel) plans as described in the previous section. The target prescription dose D_pres is 70 Gy. All the treatment plans were normalized so that 95% of the PTV received the prescription doses. 100% of the ITV and GTV received the prescription dose. All procedures are coded in the C++ programming language, and IBM ILOG CPLEX Callable library is used to solve the optimization model described above. By successively adjusting the penalty weights of all the structures associated with individual dose-volume goals in the objective function, the Anatom, FuncAvoid, and Voxel treatment plans are alternately solved to Abbreviations: DVH, dose-volume histogram; PTV, planning target volume.
improve target coverage and OAR sparing. The adjustments for each model are also guided by and made to emulate the best plan results from the other model. The process terminates when successive improvement of plan quality for either model is sufficiently small. This approach of plan generation is similar to what was used by Rao et al. [16] in an effort to make the ensuing plan comparison as fair as possible. The final dose distributions for both plans are normalized so that 95% of the PTV receive the prescription dose.
An Anatom plan is achieved by equal weighting of all the structures in the objective function. Three structures are included in FuncAvoid plan, which are high functional, medium functional, and low functional region. Functional voxels are normalized between 0 and 1. Low, medium, and high functional structures are in the region of 0 ∼ 33%, > 33% ∼ 66%, and > 66% ∼ 100%, respectively. Voxel plan is achieved by adjusting weights for lungs as a function of ventilation mapping. Weights are adjusted by trial and error methods for all the plans.
Fixed field IMRT technique was used with 6 MV Varian TrueBeam machine. This is a retrospective study. We used the exact same clinical beam setup following clinical protocols, and the exact same fields were used as the anatomic treatment plan in the clinic.

Statistical analysis
To compare the Anatom,FuncAvoid,and Voxel treatment plans, dose-volume histograms (DVHs) for the PTV and OARs and lung dose-function histograms (DFH) [17] are computed. DVH and DFH parameters collected for comparison include lung V 20 , V 30 , V 40 , mean lung dose (MLD), and their functional counterparts, mean heart dose, mean esophagus dose, maximum dose delivered to the spinal cord, and PTV uniformity index. The functional mean lung dose (fMLD) was computed for each plan as where f v and D v denote the local ventilation and delivered dose at voxel v, respectively, and L is the set of voxels in the lung. The PTV uniformity index is defined as the ratio of the dose providing 95% coverage of the PTV to the dose providing 5% coverage of the PTV. It is calculated using following formula given in Prabhakar. [18] Uniformity where Volume(Dose i ) is the volume receiving a given dose i, V presc is the volume of the prescription isoline, and PD is the prescription dose. Note that due to plan normalization, the dose providing 95% coverage of the PTV is always equal to the prescription dose. A paired sampled t-test is conducted to see whether there is significant difference at the collected DVH and DFH parameters between the two plans. A pvalue less than 0.05 ( -value) is considered statistically significant. Figure 1 shows the dose-volume histograms (DVHs) for the PTV, OARs, and lung dose-function histogram (DFH) of the Anatom, FuncAvoid, and Voxel Intensity Modulated Radiation Therapy (IMRT) plans averaged over 19 patients. The Voxel plan evidently achieved more favorable DVH and DFH for the total and functional lung than the Anatom and FuncAvoid plan. With similar PTV coverage, the PTV uniformity index is 1.07 for anatomical, and 1.08 for both FuncAvoid and Voxel treatment plans. Figure 2 shows the mean dose comparison of Anatom, FuncAvoid, and Voxel treatment planning for cord, heart, esophagus, and PTV. The dose delivered F I G U R E 1 Dose-volume histogram (DVH) for the planning target volume (PTV) and organs-at-risk, and dose-function histogram (DFH) for the lung.

F I G U R E 2
Comparison of mean doses of cord, heart, esophagus, and planning target volume (PTV) among the anatomical, functional-avoidance, and voxel-based planning.
to the heart, esophagus, and PTV are similar in all three plans. Only cord gets slightly worsened in the Voxel treatment planning compared with the other two treatment plans. We also conduct paired sample t-tests among all three treatment plans for PTV, heart, esophagus, and cord. f 2 shows the t-test results. One thing to notice is that PTV has no mean differences among all the three plans. Other structures have a mix  Lung DVH and DFH parameters in the Anatom, Fun-cAvoid, and Voxel plans for all patients are summarized in Figure 3. In this figure, the units in the Y-axis are the actual mean dose for MLD and fMLD, while it is the volume receiving the dose for V 20 , V 30 , and V 40 . Voxel plan led to an overall significant decrease in lung dose compared to Anatom and FuncAvoid planning. For total lung, the average reductions in FuncAvoid plan from Anatom plan, Voxel plan from Anatom plan, and Voxel plan from FuncAvoid plan are 7.0%, 16.8%, and 10.6%, respectively, for V 40 ; and 0.4%, 6.4%, and 6.0%, respectively, for mean lung dose (MLD). For functional lung, the reductions are 8.8%, 17.2%, and 9.2%, respectively for fV 40 ; and 1.1%, 6.2%, and 5.2%, respectively for functional mean lung dose (fMLD). These average volume reductions are shown in Table 4. We conduct paired sampled t-tests for V 20 , V 30 , V 40 , and MLD with their dose-function counterparts, fV 20 , fV 30 , fV 40 , and fMLD for all three plans. Table 3 shows the t-test results. The reductions in functional lung are statistically significant.
We also conduct paired-sample t-tests for V 20 , V 30 , V 40 , MLD, fV 20 , fV 30 , fV 40 , and fMLD to see whether there is any statistical significance among the three plans. Table 4 shows the t-test results between Fun-cAvoid with Anatom, Voxel with Anatom, and Voxel with FuncAvoid plans.Almost all the pairs are significantly dif-

DISCUSSION
This study demonstrates that the incorporation of 4DCT-based ventilation imaging in IMRT treatment planning led to significant reductions in lung dose without significantly impacting doses to other organs-at-risk (OARs) or planning target volume (PTV) coverage.
A number of studies [11][12][13][19][20][21][22][23] investigate the utilization of functional imaging in radiotherapy treatment planning (RTP). These studies generally involve different treatment techniques (e.g., 3D conformal radiation therapy and IMRT) and/or functional imaging modalities (e.g., single photon emission computed tomography (SPECT) perfusion imaging, magnetic resonance TA B L E 4 Paired-sample t-test and% of volume reduction for different volumes of total lung and functional lung among three treatment plans (MR) and 4DCT-based ventilation imaging). In addition, depending on how the regional function information is utilized, the optimization methods are rather different. For example, Seppenwoolde et al. [19] optimize the treatment plan by minimizing mean perfusion-weighted lung dose derived from the perfusion-based dose-function histogram, while Christian et al. [20] consider the planning objective to minimize the dose to a volume of functional lung contoured on the perfusion map.
In spite of these diversities, all the studies demonstrated that the incorporation of functional imaging into the RTP could improve functional sparing of the lung, at least for a particular subset of patients (e.g., patients with large perfusion deficit). The improvement was generally reflected by significant reductions in some dose-function parameters on the total lung [13,19] or certain dose-volume parameters on the functional lung). [12,20] However, there are some inconsistencies regarding dose-volume parameters on the total lung. For example, Bates et al. [23] show significantly reduced V 20 consistent with this study, while Yaremko et al. [11] reported increased total lung V 20 .
In addition, although it is generally shown by these studies that Voxel and FuncAvoid treatment of the lung could be achieved without significantly increasing doses to OARs, such as heart and esophagus, an inconsistency in the PTV dosimetric was observed. For example, Yamamoto et al. [12] report significantly degraded PTV uniformity consistent with the present study, while McGuire et al. [21] show comparable uniformity.
Esophagus dose was set in between 30 and 34 Gy, and it was 40 to 45 Gy for Cord in our optimization model. In Table 2, we see that Voxel-based plans have slightly higher cord Dmax and esophagus Dmean (p < 0.001). This means that some doses were shifted to the critical structures like the cord and esophagus. However, overall, the doses were still well below the dose constraint settings for OARs. The results of this research showed that this type of advanced treatment planning approach should be utilized in clinical settings. In clini-cal treatment planning, due to the size of the non-target non-OAR normal tissue, constraints were not imposed on these normal tissues. In this research, we added another normal tissue constraint to limit the dose shifting from lung region to the critical organs and other non-target non-OAR normal tissues. Because we are using an exact mathematical formulation, our approach guarantees the optimal solution under the current formulation. Without being able to directly handle DFH based objective and constraints, due to the limitations of the current treatment planning system, especially the solution approach behind it, we noticed sometimes the anatomical plans generated were not pareto optimal. In addition to the advantage of volume functional areas, we noticed that the overall treatment quality plan can also be improved when voxel-based optimization techniques are used. That is the additional advantage of our approach.
We now focus on the comparison of our results with those reported in St-Hilaire et al. [13] due to similar optimization approaches. They conclude functional planning to be more effective in reducing dose to the low-dose regions of the lung. However, our results show that Voxel treatment planning outperforms both Anatom and FuncAvoid plans in the case of total lung and functional lung. Table 4 shows the reduction of volumes from Anatom plan to FuncAvoid and Voxel plan and FuncAvoid to Voxel plan. The relative improvements in V i and fV i parameters increase for i = 20, .., 60. Even though improvement decreases after i = 40, the improvements in FuncAvoid over Anatom plan for V 20 and fV 20 are insignificant, which is shown in Table 4. Voxel plans reduce volume significantly in all cases compared with both Anatom and FuncAvoid plans.
From our experience we see that voxel-based plans usually use 2%-4% more MUs than anatomic plans. This is due to the advanced optimization techniques used so that the plan requires slightly more MUs to deliver a higher quality plan. We are not conducting leaf sequencing in this research; therefore, the number of segmentations is not relevant in this research.

CONCLUSIONS
This research uses a column-generation approach to generate treatment plans for lung cancer patients. Three treatment plans are generated using real patient data. Results show that Voxel treatment plan improves volume reduction compared with FuncAvoid and Anatom treatment plans.

AC K N OW L E D G M E N T S
The authors would like to thank the anonymous referees for their extremely thorough review, which led to significant changes that substantially improved the quality of this paper.

C O N F L I C T O F I N T E R E S T
The authors have declared no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available on request from the corresponding author. The data are not publicly available due to privacy or ethical restrictions.