Determination of an optimal organ set to implement deformations to support four‐dimensional dose calculations in radiation therapy planning

Surface‐based deformable image registration to generate a four‐dimensional (4D) dose calculation in radiation treatment planning requires the selection of a set of organ contours representing a basis set from which to generate anatomic deformation. The purpose of the present work was to determine the optimal set of organs needed to generate a basis set for deformation in treatment planning for thoracic tumors, such that the required computations are minimized, but that dose accuracy is high. Using retrospectively reviewed records and a deformable model algorithm in a research version of a commercial radiation treatment planning system, we calculated 4D dose distributions based on treatment plans for 10 patients with thoracic tumors. Various combinations of organs (total lungs, heart, spinal cord, external body surface) were used to generate the basis set used in the calculations for deformations. The external surface contour did not have a noticeable effect on the dose calculation. Total lung, heart, and spinal cord together provided an adequate set of deformation organs to generate accurate dose deformations. The magnitude of the calculated dose differences had no obvious relationship to tumor parameters, including site, histologic type, disease stage, extent of motion, or degree of centralization. PACS numbers: 87.55.D‐, 87.55.dk, 87.55.kh


I. INTRODUCTION
The anatomic deformation caused by respiratory activity constitutes an obstacle to the accuracy of dose calculations for radiation therapy of tumors located in the thorax. By explicitly imaging tumor movement, four-dimensional (4D) computed tomography (CT), which involves the generation of sequential three-dimensional (3D) image data sets representing phases of the patient's respiratory cycle, has provided greater accuracy in the radiation treatment of thoracic tumors. (1) However, techniques that would account for motion-induced deformation of tumor and surrounding tissues have not been well studied or clarified. Accurate calculation of the radiation dose distributions in a 4D CT data set requires a determination of the manner in which the 3D CT image of a reference phase deforms to all other phases in the 4D data set.
One technique for determining the deformation of a 3D CT data set is to track the deformation of specified anatomic structures from phase to phase in the 4D data set. (2,3) In this approach, a tetrahedral mesh representing the surface of an anatomic structure is generated from contours a drawn in the radiation treatment plan, and basis functions are used to interpolate deformations from the surface components through the volume of the known structure. This technique was recently incorporated into a research version of a commercial radiation treatment planning system (Pinnacle 3 , v8.1t: Philips Medical Systems, Milpitas, CA). (4) The current implementation of 4D dose calculation requires that the user identify a set of organs that will undergo deformation. Points other than those explicitly identified in the mesh are interpolated using an elastic body spline deformation model, which characterizes the deformation of elastic bodies under the influence of elastically applied forces. (4) The treatment planning system gives the user some flexibility in the selection of the organ surfaces that will represent the deformation of the thoracic cavity. In making that selection, the user must strike a balance between accuracy and efficiency. That is, the greater the number of organ surfaces selected, the more accurate the dose calculation is likely to be; however, because the chosen organs must be delineated onto the 3D data sets and mathematically deformed to simulate phases of respiration (a time-consuming process), the minimum number of organ contours necessary for dose accuracy should be used in the calculations. The full set of deforming anatomic structures that can be used for this purpose consists of the lungs, heart, spinal cord, and external body surface. However, it may be that a subset of these organ contours can be utilized in dose calculations without unduly compromising dose accuracy. In the current study, we examined that possibility and also looked for correlations between tumor characteristics (size, histologic type and disease stage, degree of centralization within thoracic cavity, and extent of motion) and sensitivity to the organ contours used to approximate the thoracic deformation.

A. Patient selection
In this institutional review board-approved study, we retrospectively reviewed the records of 10 patients who had undergone radiation treatment for tumors in the lung (n = 8) or esophagus (n = 2) and for whom 4D CT data sets had been acquired. Patient treatment plans were selected sequentially from the patient database in the treatment planning system and covered the period June and July 2006. Among the lung tumor patients, 1 received hypofractionated radiation treatments; the remaining 7 lung tumor patients and both patients with esophageal cancers received radiation using conventional fractionation schemes. Table 1 summarizes tumor site, disease stage, and prescription dose for the 10 patients.

B. Treatment plans
Patients were immobilized using a conventional vacuum-bag immobilization technique and were scanned using a multislice helical CT scanner (Discovery ST: GE Healthcare Systems, Waukesha, WI) operating in cine mode. The 4D CT image data sets, consisting of 3D CT data sets for each of 10 equally-spaced phases of the respiratory cycle, were obtained using an image-binning approach. (1) Respiration was monitored using an external marker resting on the patient's abdomen (RPM: Varian Medical Systems, Palo Alto, CA). The CT images were reconstructed with a slice thickness of 2.5 mm. We planned the treatments using the treatment planning system mentioned earlier, determining an explicit internal target volume (ITV) by mapping the motion of the clinical target volume (CTV) on the 4D CT data set. Tumor targeting was based on information in the 4D CT data set and was transferred to either a free-breathing CT data set or an average CT data set, on which the treatment plan was developed. (An "average data set" is defined as a 3D CT data set in which each voxel value is the average value among the 10 phases that make up the 4D CT data set.) In generating the planning target volume (PTV) for conventional fractionation plans, we established an isotropic 1.0-cm margin around the ITV to account for setup uncertainty and interfractional variability in respiration. Conventionally fractionated patients were typically planned using inverse planning for intensity-modulated radiation delivery with 5 -7 6-MV photon beams. For image-guided hypofractionated treatments, the margin used to generate the PTV was 0.3 cm around the ITV. Patient immobilization was achieved with an extended vacuum bag, and image guidance was provided by daily in-room CT imaging immediately before irradiation. Planning was accomplished using fixed 6-MV photon beams. Anatomic structures that were normally contoured included the spinal cord, the right and left lungs, total lung, esophagus, and heart. Beam configurations were determined in a variety of ways, including conventional 3D conformal radiation therapy, intensity-modulated radiation therapy, and image-guided hypofractionated radiation therapy.
The treatment plans-including all CT data sets (planning set plus 10 phases comprising the 4D data set), the contours of the PTV and normal anatomic structures, and the beam configurations-were copied into a research directory, and the patients' identities were masked before further analysis was done.

C. 4D dose calculations
To mathematically represent the deformation of the thoracic cavity throughout the respiratory cycle, we used the surface-based deformable image registration method. (2,3) Regions of interest (ROIs), which may be organs or target volumes, were delineated on a single reference 3D CT image data set either by the attending radiation oncologist or by the treatment planner. (In the latter case, plans were subsequently reviewed by the radiation oncologist.) A deformable model algorithm was then used to propagate the ROIs to the phases constituting the 4D CT image data set. (5,6) The 4D dose distributions were calculated first by using the collapsed cone convolution algorithm to perform a dose calculation on each of the 10 phases of the 4D CT data set. (7) One data set, typically the 0% phase data set, was identified as the reference data set, and each of the other data sets was deformed to the reference data set using the deformable image registration technique (4) based on a specified set of ROIs and the elastic body spline deformation model. The dose matrix in each data set was then deformed to the reference data set based on the foregoing deformations, and the doses were calculated and accumulated.
Initially, the ROIs used for the deformation included the total lung (L), heart (H), spinal cord (C), and external body surface (E). This LHCE set of ROIs represented the highest attainable dose accuracy for the 4D dose calculation, because it included all anatomic structures that are routinely delineated in the thoracic region for treatment planning. We then recalculated the deformation using subsets of the control organ set by excluding selected contours from the deformation. The following sets of ROIs were used in the analysis: LHCE (control), LHC, LH, LC, LE, and L. In all cases, the total lung was included because it represents the major contribution to deformation.

D. Dose-volume histogram analysis
First, we used the LHCE set of ROIs to generate control dose-volume histograms (DVHs) for each patient. Then, we used each of the subsets of ROIs to generate experimental DVHs based on 4D dose calculations. Finally, we compared the experimental DVHs with the control DVHs. We generated DVHs for the spinal cord, CTV, PTV, heart, right lung, left lung, and both lungs. Several parameters were selected to analyze the DVHs, including the V20 value (fractional volume of the ROI receiving at least 20 Gy) for the total lung ROI, the D5 value (dose received by 5% of the ROI volume) for the spinal cord, the D95 value (dose received by 95% of the ROI volume) for the PTV, and the V(prescription) value (fractional volume of the ROI receiving the prescription dose) for the CTV.
Total lung V20 was selected because it appears to be closely related to the risk of radiation pneumonitis. (8) Although the maximum dose to the spinal cord is used to predict the risk of radiation myelitis, (9) we studied the cord D5 value instead because the maximum dose to the cord may vary greatly for situations in which the DVH exhibits a long "tail." (The D5 value is likely to be more consistent among deformation organ sets.) We selected D95 for the PTV because our clinical standard of practice is that 95% of the PTV receives the prescription dose. Finally, we selected V(prescription) for the CTV because our clinical standard of practice is that the entire CTV receives the prescription dose [that is, that Deviations of the experimental DVHs from the control DVHs were visually characterized by means of composite DVHs and were then quantified using the differences between the D95 values of each organ subset and the control (LHCE) organ set. We organized the results by extent of deviation from control and then looked for trends that might be related to tumor size, histologic type, or disease stage; approximate degree of centralization; and extent of respiratory motion, which was approximated by measuring the range of motion of the apex of the right side of the diaphragm.

III. RESULTS
As a general rule, the DVHs for total lungs were similar, regardless of the organs used in the deformation. Fig. 1(a) shows the DVHs for total lung for patient 5, who showed the largest difference in DVHs calculated using various organ sets; Fig. 1(b) shows the DVHs for total lung for patient 7, whose results were more typical of the entire group of 10 patients. Both figures demonstrate that the DVHs for each ROI set were quite similar, with the exception of the DVH based solely on deformation of the lung contours. Table 2 summarizes the V20 values, showing deformations based on each combination of ROIs. These data indicate that the LHC deformation organ set yielded V20 values that most closely approximated those obtained with the full LHCE set [root mean square (RMS) deviation: 1.25%; maximum deviation: 2.59%], although the LH and LE deformation sets yielded comparable values. In general, deformation based on L alone gave the poorest agreement with LHCE, although the LC sometimes yielded even more divergent results. Table 3 shows the data from Table 2 rearranged in descending order of absolute difference between the LHC-and LHCE-generated total lung V20 values. Tumor site, disease stage, prescription dose, and hypofractionation status are also indicated. Visual inspection of the data detected no trends linking the differences in V20 values with tumor site, disease stage, or fractionation scheme.  Table 4 shows values of D5 and differences between D5 values calculated using each subset of deformation organs and the control set. Again, calculations based on the LHC subset had the smallest RMS deviation (235 cGy) and the smallest maximum deviation (500 cGy). Fig. 2(a) shows the cord DVHs for patient 5, for whom the disagreement between calculations was large, and Fig. 2(b) shows the cord DVHs for patient 7, whose results were more typical of the entire patient group. The cord dose calculated for patient 5 using the LHC deformation organ set would most likely be viewed as clinically acceptable, because a negligible amount of cord is calculated to receive a dose greater than 45 Gy, but the cord doses calculated using the LH and LE data sets likely would not be considered clinically acceptable because a non-negligible amount of cord is calculated to receive a dose greater than 50 Gy. TABLE 3. Patients sorted in descending order by absolute value of the difference between total lung V20 obtained by using lung, heart, cord (LHC) for deformation and total lung V20 obtained by using lung, heart, cord, external contour (LHCE) for lung deformation   Table 5 shows the data in Table 4 rearranged in descending order of the absolute difference between the LHC-and LHCE-generated cord D5 value, together with tumor site, disease stage, prescription dose, and hypofractionation status. Although the absolute D5 differences for esophagus treatments are clustered, the sample is too small for that clustering to be considered any sort of a trend. The difference in D5 values for the hypofractionated treatment was 0, but that result may simply have been the result of the spinal cord receiving very little dose in the hypofractionated treatment. Table 6 compares values of D95 for the PTV, together with the differences in D95 between the control calculations and those derived using each subset of deformation organs. Once again, the LHC subset exhibited the smallest RMS deviation (232 cGy) and the smallest maximum deviation (500 cGy). For many patients, the differences in the dose calculated using only the lungs as the deformation organ set and the dose calculated using the full LHCE organ set were so great that the L subset dose distribution would be viewed as clinically unacceptable.
Finally, Table 7 shows the values of V(prescription) for the CTV and the differences in V(prescription) between calculations performed using each subset of deformation organs and those performed using the LHCE set of organs. In some cases, the volume receiving the actual prescription dose was somewhat smaller than originally intended, and thus the volume receiving the lower dose is recorded here. Again, the LHC subset exhibited the smallest RMS deviation (4.59%) and the smallest maximum deviation (11.88%). In several cases, use of the total lung as the only ROI for deformation resulted in severe underdosing of the CTV. Figs. 3(a) and 3(b), for patients 5 and 7 respectively, vividly illustrate that result. The data sets for patient 5, which exhibited the largest disagreement between other organs, also had a large (but not the largest) disagreement for the CTV; the degree of the disagreements between the data sets for patient 7 were more typical of those for the whole group. For patient 5, the only organ subset that gave a DVH for the CTV comparable to that calculated using the LHCE organ set was the LHC set; DVHs calculated using the other sets suggested that the doses were clinically unacceptable. For patient 7, DVHs calculated using the LHC and LE organ sets were comparable to those from the LHCE calculations, and DVHs derived using LC and LH indicated doses that likely would be considered clinically acceptable; however, DVHs calculated using L alone were clinically unacceptable.

IV. DISCUSSION
To calculate 4D dose distributions, a Cartesian dose calculation matrix is generated on each phase in the 4D CT data set, the dose distribution on each of the dose matrices is calculated, and then each dose matrix is deformed to the reference phase, accumulating the dose on the reference matrix for that phase. The difficult part is the deformation of each phase, because the matrices that define the deformation ideally indicate how each tissue volume element moves between phases during respiration. Current approaches to deformable registration fall into two classes: image-based deformations and model-based deformations.
Image-based deformations, as exemplified by algorithms such as the "demons" algorithm (10) and optical flow methods, (11) generate deformations based on the values of the voxels in the CT image data sets. Although user intervention may not be needed to effect those deformations, one-to-one correspondence between CT voxels and tissue volume elements is not guaranteed. In an extreme case, in which deformation of a homogeneous medium occurs, image-based deformable registration methods may not be able to detect the magnitude or even the presence of a deformation.
Model-based deformations, such as that described by Kaus et al. (4) and implemented in the commercial radiation treatment planning system used here, require pre-processing in the form of segmentation of each data set to generate the contours of corresponding ROIs. Additional information in the form of basis functions and biomechanical information (12) may have to be introduced to interpolate deformations into the interior of segmented anatomic structures. In addition, some judgment is needed to select a set of organ contours sufficient to accurately reflect the totality of deformations that the interior of a patient undergoes during the respiratory cycle.
On the basis of our findings, we conclude that accurate doses can be calculated in 4D CT planning of radiation treatments using primarily internal organ contours for image deformation, rather than the full set of contours typically used, including those for external body surface. As seen in the present work, selecting simply the total lungs as the deformation organ does not adequately characterize the deformation needed to accurately calculate the 4D dose distribution. It appears that organs such as the heart and spinal cord, medial structures that do not deform as extensively as the lungs, are needed to calculate dose accurately. However, including more organs than is necessary slows down the deformation process in two ways. First, organ contours have to be generated on each phase of the 4D data set. That task is accomplished by copying the mesh model of the organ from the reference CT data set to each of the other phases, and then by deforming the model to fit the specific images on the phases. Second, the deformations of all organs have to be incorporated into the deformation calculations. Consequently, a good compromise needs to be found to limit the number of organs included in the deformations.
In the present work, we found little difference between DVHs generated using the original organ set (LHCE) and those generated using the experimental subset LHC, suggesting that the external body surface contour does not necessarily play a large role in characterizing data set deformations. Calculations based on deformations using other organ subsets tended to exhibit substantial variation, with the DVH generated using the lung contour alone generally being the least accurate.
Tables 3, 5, and 8 reveal no obvious trend linking the differences in the DVHs generated using the LHC organ subset and the full LHCE set with the tumor site or disease stage. Patient 6, who received hypofractionated radiation treatment, had a larger-than-typical discrepancy between CTV dose calculations; however, that patient had less discrepancy in dose calculations for other ROIs. For most other patients, there appeared to be correlation between DVH consistencies among ROIs.
Reducing the number of organs used for a deformation does not always result in a less accurate calculation. For instance, the doses calculated for Patient 1 appeared to depend solely on the deformation of the lungs. The reason for this complete insensitivity of dose to deformation of the heart, spinal cord, and external contours is unclear.
We further found in a partial study that inclusion of the GTV as a deforming organ does not appear to affect the dose calculation.
The substantial increase in the accuracy of predicted radiation doses to thoracic tumors afforded by 4D dose calculations may result in more effective treatment planning. As noted here, a major limiting factor in the practicality of the 4D approach, the volume of data calculations required, may be reduced if further studies can determine precisely which tumor parameters affect the sensitivity of the dose calculation to the organ subset. If tumor characteristics can be correlated with the organ subsets required for an accurate dose calculation, then it may be possible, through the visualization of CT scans alone, to individually tailor the organ sets used in each treatment plan to generate a dose calculation with an ideal balance between accuracy and efficiency.
Finally, it should be noted that the present study has demonstrated that the LHC subset of organs gives the same dose distributions as the complete LHCE set of organs, but does not make any statements regarding the absolute accuracy of 4D dose calculations. Comparison with measurements would verify the accuracy of 4D dose calculations, and such work is in progress.

V. CONCLUSIONS
The present work has demonstrated that, for 4D dose calculations in the thorax using a surfacebased deformable image registration approach based on deforming contours of specific organs, a minimum set of total lung, heart, and spinal cord is necessary to give results comparable to calculations using all organs.

ACKNOWLEDGMENTS
The present work was supported in part by a Sponsored Research Agreement from Philips Medical Systems. The authors also acknowledge the assistance of Dr. Veni Ezhil and Michael Kantor.