Synthetic 4DCT(MRI) lung phantom generation for 4D radiotherapy and image guidance investigations.

PURPOSE: Respiratory motion is one of the major challenges in radiotherapy. In this work, a comprehensive and clinically plausible set of 4D numerical phantoms, together with their corresponding "ground truths," have been developed and validated for 4D radiotherapy applications.
METHODS: The phantoms are based on CTs providing density information and motion from multi-breathing-cycle 4D Magnetic Resonance imagings (MRIs). Deformable image registration (DIR) has been utilized to extract motion fields from 4DMRIs and to establish inter-subject correspondence by registering binary lung masks between Computer Tomography (CT) and MRI. The established correspondence is then used to warp the CT according to the 4DMRI motion. The resulting synthetic 4DCTs are called 4DCT(MRI)s. Validation of the 4DCT(MRI) workflow was conducted by directly comparing conventional 4DCTs to derived synthetic 4D images using the motion of the 4DCTs themselves (referred to as 4DCT(CT)s). Digitally reconstructed radiographs (DRRs) as well as 4D pencil beam scanned (PBS) proton dose calculations were used for validation.
RESULTS: Based on the CT image appearance of 13 lung cancer patients and deformable motion of five volunteer 4DMRIs, synthetic 4DCT(MRI)s with a total of 871 different breathing cycles have been generated. The 4DCT(MRI)s exhibit an average superior-inferior tumor motion amplitude of 7 ± 5 mm (min: 0.5 mm, max: 22.7 mm). The relative change of the DRR image intensities of the conventional 4DCTs and the corresponding synthetic 4DCT(CT)s inside the body is smaller than 5% for at least 81% of the pixels for all studied cases. Comparison of 4D dose distributions calculated on 4DCTs and the synthetic 4DCT(CT)s using the same motion achieved similar dose distributions with an average 2%/2 mm gamma pass rate of 90.8% (min: 77.8%, max: 97.2%).
CONCLUSION: We developed a series of numerical 4D lung phantoms based on real imaging and motion data, which give realistic representations of both anatomy and motion scenarios and the accessible "ground truth" deformation vector fields of each 4DCT(MRI). The open-source code and motion data allow foreseen users to generate further 4D data by themselves. These numeric 4D phantoms can be used for the development of new 4D treatment strategies, 4D dose calculations, DIR algorithm validations, as well as simulations of motion mitigation and different online image guidance techniques for both proton and photon radiation therapy.

INTRODUCTION

Anatomical motions are a challenge for many medical applications. In particular, respiratory motion is a major issue in radiation therapy, as it can negatively and substantially affect the quality of treatment. For imaging such motions, 4DCT is currently considered as the clinical standard. , However, retrospectively reconstructed 4DCTs only represent an average breathing cycle and cannot capture irregular breathing motion. Indeed, irregular breathing during acquisition can lead to substantial artifacts in the resulting 4D data sets. Alternatively, 4DMRI can be acquired over long periods and can help to capture the irregularity of the breathing pattern. This modality, however, does not provide information on tissue density, which is essential for many radiotherapy and other imaging applications. Anatomical motion is also a complex, multivariate problem that is not always reproducible and is often unpredictable, making the ability of 4D imaging techniques to capture the true motion questionable. Much has been reported in the literature on imaging and mitigating motion. However, due to the complexity of the problem, it is often extremely challenging to validate such techniques, as in real patient data, there is no ground truth knowledge of the true, underlying motion that is being imaged or mitigated. For this reason, numerical 4D phantoms have been proposed and widely used in research. The motion of a numerical phantom, even though it might not exactly reflect the motion of real patients, is known exactly and can therefore be considered as the “ground truth” motion of the phantom and used for algorithm validation.

Many such numerical phantoms have been developed. A detailed description of existing phantoms and their applications can be found in comprehensive review articles. , In very early examples, organs are represented analytically as simple geometrical forms like spheres, ellipses, cylinders, or planes. One such example is the Medical Internal Radiation Dose (MIRD) phantom representing an average adult male. Such stylized phantoms can provide a lot of flexibility to model organ motion or variable geometries but are not detailed or realistic enough for assessing complex processes like 4D imaging or treatment planning in radiotherapy. On the other hand, voxelized phantoms based on segmented image data can give more detailed representations of anatomy, but their flexibility is limited in terms of motion and geometric variability. Two examples are the early VIP‐man or the Gesellschaft für Strahlenforschung (GSF) phantom. Later, hybrid phantoms, combining the complexity and realism of voxelized phantoms with the flexibility of stylized phantoms, have also been proposed, for example, the XCAT phantoms , , , and the 4D VIP‐man. These phantoms are based on patient imaging data (CT, MRI), and the surfaces of organs are described by polygonal meshes or with the help of non‐uniform rational B‐splines (NURBS). In the original XCAT phantoms, each organ was modeled as being uniform, apart from some modeled airways in the lung, but some work has also been put into modeling more detailed and realistic structures in the lung and bones. , The XCAT phantoms also further include cardiac and respiratory motion models. Respiratory motion is controlled through a curve describing the diaphragm motion in superior–inferior (SI) direction and one controlling the anterior–posterior (AP) expansion of the chest.

As they have developed, numerical motion phantoms have been increasingly used in radiation therapy to assess and validate 4D imaging, online image guidance and motion mitigation techniques. In the context of 4D imaging, phantoms were used for the investigation of image reconstruction algorithms or mitigation of motion artifacts for 4DCT, 4D cone‐beam CT (CBCT), , , 4D positron emission tomography (PET), , , and single‐photon emission computed tomography (SPECT). Likewise, phantoms have been used to simulate online image guidance scenarios, such as digitally reconstructed radiographs (DRRs), electronic portal imaging, surface motion, cine MRI, or ultrasound guidance. Phantoms also provide the “ground truth” motion and can therefore be used for motion modeling , , as well as for deformable image registration (DIR) performance validation. In photon radiotherapy, phantoms have found applications in treatment planning studies and simulations of gating and tracking. , , ,

Arguably, however, the most sensitive radiation therapy modality to motion is pencil beam scanned (PBS) proton therapy. Due to the sharp distal fall‐off of the proton depth–dose curve and the sequential (in time) application of individually weighted and narrow proton “pencil beams” (each just a few millimeters in diameter), PBS proton therapy is especially sensitive to patient/tumor geometry changes and the interference with tumor motion. , , As such, bespoke numerical phantoms have been used to study various motion mitigation strategies in proton therapy, such as rescanning, gating, , or tracking or for the development of new 4D dose calculations or 4D planning strategies. As such, the wider availability of flexible and anthropomorphically accurate numerical phantoms, for which “ground truth” motions are known, and which also accurately represent both realistic internal motions and densities, would be of great value to both the radiotherapy and imaging research communities as benchmark data sets for assessing motion imaging and mitigation techniques.

In this paper, we describe a comprehensive set of numerical 4D lung phantoms that have been developed at our institute. For the phantom generation, we expand on an approach already reported by our group, , which mapped motion from 4DMRIs onto static reference CTs using manual landmarks to establish the correspondence between the CT and a reference phase of the 4DMRI. More recently, a similar approach has been reported using multimodality DIR to establish such correlations. In this present work, however, we apply DIR to the segmented lung masks of each subject for high quality inter‐subject geometric correspondence, enabling patient geometries (CTs) to be animated with a large variety of realistic deformable motion patterns (from 4DMRIs) to create an extensive 4D data set for radiotherapy and imaging applications. These patterns also provide the “ground truth” motions for the data set. All developed data sets have been validated using 4D dose calculations for PBS proton therapy which, as described above, we believe is one of the most sensitive radiotherapy modalities to both motion and density variations.

To begin with, we first introduce the image data set preparation procedures and the workflow used to develop the set of 4D numerical phantoms. After the geometrical validation of the phantoms and the evaluation of the densities, an assessment of the accuracy of the motion and modeled density changes within the phantoms is performed using 4D PBS proton therapy dose calculations.

MATERIALS AND METHODS

The numerical lung phantoms developed in this work are based on real patient geometries from CTs, combined with motion information from 4DMRIs in order to create synthetic 4DCTs, also called 4DCT(MRI)s. We applied the proposed workflow to 13 different CTs from lung cancer patients, each with pronounced differences in lung shapes/volumes and tumor sizes/locations. Motions have then been extracted from multi‐breathing‐cycle 4DMRIs of five volunteers using DIR. All motion data used for the generation of the phantoms (in the form of moving lung meshes) can be accessed via the following link: https://doi.org/10.5281/zenodo.5016294, while the source code for generating 4DCT(MRI)s can be downloaded from https://doi.org/10.5281/zenodo.5010964.

Base‐line data sets

Clinical lung cancer CT data

CTs of 13 different lung cancer patients have been used to provide a wide range of realistic anatomical geometries for our phantoms. CT1–CT6 are end exhale (EE) 4DCT phases of six early‐stage non‐small‐cell lung cancer (NSCLC) patients, whereas CT7–CT13 were acquired using visually guided voluntary deep‐inspiration breath‐hold (DIBH) for seven advanced‐stage NSCLC patients. , All CTs were acquired on Siemens CTs (Siemens Healthineers, Erlangen, Germany) with a resolution of 0.98 × 0.98 × 2 mm3. The 13 CTs thus exhibit large variabilities in tumor size, tumor location, lung shape, and volume. For all CTs, both lungs were segmented manually by physicists. The 3D rendering of the lungs and corresponding volumes are shown in Figure 1.

FIGURE 1

(a) Top: CT geometries from 13 patients. CT1–CT6 are end exhale (EE) reference CTs extracted from 4DCTs, whereas CT7–CT13 were acquired during deep‐inspiration breath‐hold (DIBH). Bottom: MRI lung geometries from five volunteers at EE reference phases from the 4DMRIs. (b) Volumes of both lung halves for CTs and MRIs. For the 4DMRIs the volumes for both the EE and the end inhalation (EI) reference state are shown

Volunteer lung 4DMRI data

The time‐resolved deformable respiratory motion was derived from time‐resolved 3D volumetric MRI data (4DMRIs) of five volunteers acquired on a 1.5 T MRI scanner (MAGNETOM Aera, Siemens Healthineers). The volumetric MRI had a 500 × 275 × 400 mm3 field of view with a 3.125 × 3.125 × 3.125 mm3 image resolution. For imaging, a motion‐aware acquisition technique was used which is based on a spoiled gradient echo sequence which acquires core and pseudo‐randomly sampled peripheral k‐space patches alternately. The motion‐aware reconstruction “corrects” the spatial offset between a peripheral patch and a reference EE core utilizing DIR. This yields an effective time resolution of 2.25 Hz. MRI data were acquired over 11 min. The 4DMRI data sets used in the present study were acquired and processed as part of our previous work, to which the interested reader is referred to for further details and analysis. Both lungs were manually segmented on the reference EE and the following end inhalation (EI) phase of the reconstructed 4DMRI (EE lung shapes visualized at the bottom of Figure 1). For the generation of the example 4D phantoms, five breathing cycles following the reference EE phase for each 4DMRI were considered. In addition, a longer pattern of 20 cycles was extracted for MRI4, as an example representing breathing irregularities over a longer time‐scale. More cycles from MRI1, MRI2, and MRI5 exhibiting large motion irregularities were also selected. In total, eight deformable motion patterns consisting of 67 breathing cycles were selected.

4DCT(MRI) workflow

The originating CT data sets consist either of EE reference phases from 4DCTs or static DIBH CTs. To correctly map the motion from the 4DMRI onto the reference CT, a similar 4DMRI reference state first needs to be defined. For the EE CTs, the EE states used for the 4DMRI reconstruction and DIR were selected. For DIBH CTs, the next EI state after this 4DMRI EE state was used. A schematic of the complete workflow is illustrated in Figure 2. First, both lungs were manually segmented on both the reference CT and the reference MRI. From the segmented MRI lungs, surface meshes were extracted using the visualization toolkit (vtk).a , Second, these binary CT lung masks were then registered to the MRI lung masks using Plastimatch,b to establish anatomical correlation between the CT and MRI lung shapes. These registrations of the binary masks were performed separately for both lung halves. For all binary mask registrations, we first utilized a 3D affine registration, followed by a 3D multi‐resolution B‐spline registration using the mean‐squared error metric. The results of the lung DIR were used to deform the MRI surface meshes to match the CT lungs, such that mesh point correlation between the CT and MRI lungs is established. If the registration and deformation of the surface mesh led to unsatisfactory results (e.g., part of lung missing), strategically placed landmarks can be manually defined and used to guide another registration. These landmarks can be placed in areas where the previous registration had difficulties by selecting landmarks in similar anatomical locations on the lung surface of the CT and MRI. Moreover, grid points within the CT and MRI lungs were also inserted. For that purpose, a regular grid of points with a spacing of 15 mm was first inserted inside the CT surface mesh, which was then interpolated to positions inside the MRI mesh. Third, the Deformable Vector Fields (DVFs) resulting from registering each 4DMRI state to the reference MRI state were used to warp the reference CT mesh (including internal grid points) for generating time‐resolved CT meshes. Next, for each state the position of every CT mesh or grid point P was subtracted from the corresponding point position P in the reference mesh, and the resulting vectors (P – P starting at P ) were interpolated from the sparse deformation information to a continuous DVF for the CT using a B‐spline insight toolkit (itk)c implementation. In the same step, the DVF is also extrapolated into the abdominal region, such that the motion at the inferior border of the CT is zero unless the moving lungs would extend outside of the CT field of view, in which case the full extent of the lung motion is allowed. To preserve sliding organ motions along boundaries, the motion of the ribcage and body was set to zero by applying a mask to anatomical regions assumed to be static. On the reference CT, this was manually contoured to separate the static ribcage from the moving parts inside the thoracic and abdominal cavity. Finally, the synthetic 4DCT(MRI) data were derived by warping the reference CT with the resulting continuous DVFs.

FIGURE 2

Workflow to generate a 4DCT(MRI) based on a reference CT and a 4DMRI

It is important to note that these resulting DVFs provide the “ground truth” motion of each synthetically generated 4DCT(MRI) data set. Throughout the rest of this paper, we will use the following notation to refer to such data sets: for example CT1(MRI1) refers to the 4DCT(MRI) generated using the CT1 reference geometry and motion from MRI1.

Phantom validation

Geometrical validation

The quality of the DIR of the binary lung masks can be assessed by checking the registration results visually and quantitatively by calculating the dice coefficient. CT surface meshes were also directly generated from the segmented CT lung masks using vtk and used to visually check the similarity to the deformed CT lung meshes. The mesh point correspondence between the deformed CT meshes and the MRI meshes was visually checked in ParaViewd to make sure corresponding points are in a similar anatomical location. Further, the deformed mesh was converted into a binary mask and the dice coefficient of this deformed mask and the originally segmented CT mask for all registrations was calculated.

Density validation

For all further validations, we focused on the six 4DCTs (CT1–CT6), as they provide both density and motion information which are necessary for validation. To validate our workflow, we generated what we call a 4DCT(CT) for each of the 4DCTs by warping the reference EE phases according to the 4DCT motion itself, following the exact procedures in Figure 2. However, as the motion and the reference image came from the same subject, no DIR of binary lung masks is required. Moreover, each 4DCT phase was registered to the reference phase using a 3D multi‐resolution B‐spline registration in Plastimatch.

Based on the current procedures, density changes within the lung region are not considered in the generation of a 4DCT(MRI) or 4DCT(CT), as the positions of voxels are only warped, without changing their Hounsfield values in either the image registration or warping process. However, in reality, respiratory motion leads to a change of lung volume, whereas the mass of the lung stays approximately the same, leading to a change of the lung density. Therefore, the effect of density changes was evaluated by calculating the mean Hounsfield units (HU) in the lung (including tumor) for each phase of both the original 4DCT and the 4DCT(CT).

Additionally, we have calculated a DRR simulating AP fluoroscopic imaging for each of the six 4DCTs and the corresponding 4DCT(CT)s. Such a DRR was generated for the EI state of each 4DCT, as the biggest differences between the 4DCT and 4DCT(CT) are expected for this state due to the choice of the EE reference phase for the 4DCT(CT). The DRRs of the 4DCT and the derived 4DCT(CT) were then compared visually and in terms of relative gray level differences.

Motion and density validation

Another limitation of the current numerical phantoms is that the rib motion is not modeled. This is a consequence of masking out the body and ribs during the DIR process (to ensure smooth sliding boundaries between the lung and ribs/chest wall). To separate the effect of the ribcage motion and the density changes, we introduced another synthetic 4D image for validation, a “patched 4DCT” (p4DCT). This is obtained by combining the static ribcage of the reference phase with the changing interior regions (lungs, heart, etc.) of the original 4DCT. The static mask of the EE 4DCT reference phase, which is used for the 4DCT(CT) generation, can be applied for this purpose. For voxels in the non‐moving part the EE HU value is kept for all phases, whereas for voxels within the moving part the HU value of that voxel from the respective 4DCT phase is utilized. Therefore, the p4DCT has the same HUs within the lung as the originating 4DCT reflecting breathing‐induced density changes, but still has a static ribcage as in the 4DCT(CT).

To investigate the magnitude of the impact of these two limitations, both effects were studied using 4D dose calculations for PBS proton therapy treatment plans. Four‐dimensional dose distributions , were calculated for 3D (static), two‐field single‐field uniform dose (SFUD) plans, optimized on the planning target volume (PTV) (clinical target volume (CTV) extended by a 2 mm isotropic margin) of the reference phase of the 4DCT. The beam model and delivery dynamics of the PSI gantry 2 were used and 4D dose calculations performed using the deforming dose grid approach, which considers the different motion scenarios, as well as induced density changes from 4DCT, 4DCT(CT) and the p4DCT, respectively. For example, by comparing 4D dose distributions calculated on the 4DCT to 4D dose distributions on the p4DCT, the effect of the ribcage motion can be assessed. On the other hand, a comparison of 4D dose distributions for the p4DCT and the 4DCT(CT) indicates the effect of the density changes.

We have created static plans on the EE reference phases of CT1–CT6 with two fields. A field from the side (left–right or right–left) and a second field with a 45° angle from the front were chosen for all plans. The 4D dose distributions calculated on the 4DCT, p4DCT and the 4DCT(CT) were compared in terms of dose differences and a 3D gamma analysis using a 2% dose difference (absolute) and a distance to agreement of 2 mm (2%/2 mm). A dose threshold of 10% of the prescribed dose (2 GyRBE per fraction) was applied to define the analysis region. V 95% of the CTV and the mean dose to the healthy lung half surrounding the tumor was examined.

RESULTS

The 4DCT(MRI) numerical phantom library

The 4DCT(MRI) phantoms are based on 13 CTs and eight deformable motion patterns from five 4DMRIs. An illustration of the motion patterns is displayed in Figure 3 for one example point in the dome of each lung half (see Figure 3a). Figure 3b,c illustrates the SI amplitudes and the periods of the 4DMRI motion patterns. A total of 67 individual breathing cycles from five 4DMRIs are included, thus leading to a data set of 871 different 4DCT(MRI) breathing cycles. The animation of the CT anatomies with different motion patterns results in an average tumor center SI motion amplitude of 7 ± 5 mm. For all 4DCT(MRIS)s, averaged over the breathing cycles of each motion pattern, the minimum observed tumor amplitude is 0.5 mm and tumor displacements up to 22.7 mm were observed.

FIGURE 3

For each of the 4DMRI geometries, a point in the dome of each half of the lung was selected (see (a) for MRI1) and the superior–inferior (SI) motion is analyzed (left: red, right: green). Boxplots of amplitudes (SI) (b) and periods (c) over all breathing cycles of the 4DMRIs. SI displacements over time for MRI1–MRI5 (d–h). Selected motion patterns for 4DCT(MRI) generation are marked by gray shaded areas

Figure 4 shows example 4DCT(MRI) slices from CT2(MRI1) and CT11(MRI2). In Figure 4a–d, the EE reference phase of CT2, as well as an EI phase of CT2(MRI1), are presented. Figure 4e–h depicts the DIBH reference CT11 together with an EE slice of CT11(MRI2). Animations of one example cycle of the 4DCT(MRI)s, CT2(MRI1), and CT11(MRI2) can be found in Figure S1. Figure S2 also shows the CT2 geometry animated by a breathing cycle from the other MRIs. The variability of the CT data set is illustrated in Figure S3, where each CT is animated with the same breathing cycle following the reference phase of MRI1.

FIGURE 4

Example coronal and sagittal 4DCT(MRI) slices at different breathing phases of CT2(MRI1) (a–d) and CT11(MRI2) (e–h). For CT2(MR1) the reference end exhale (EE) phase and an end inhalation (EI) phase are shown. For CT11(MRI2) the deep‐inspiration breath‐hold (DIBH) reference and an EE phase are displayed. The clinical target volume (CTV) is contoured on each image. The white lines have been inserted for visual reference

The results of the registration of the binary masks and the subsequent deforming of the MRI lung meshes were checked visually by comparing the original CT lung mesh and the deformed MRI lung mesh. Additionally, correspondence was also visually checked for some mesh points. The dice coefficient was calculated between the originally segmented lung and the deformed lung. We obtained a mean dice coefficient of 0.963 ± 0.010 and the minimal dice coefficient is 0.931 over all investigated cases.

The mean HU values inside the lung for each breathing phase were calculated for all 4DCTs and 4DCT(CT)s. Whereas the mean HU inside the lung remains almost constant over the whole breathing cycle for all 4DCT(CT)s, a decrease of the mean lung HU between full exhalation and full inhalation can be observed for the original 4DCTs. The density changes between opposite extreme breathing phases lead to HU changes up to 49 HU for the examined 4DCTs. An average mean lung density decrease from the reference EE phase to full inhalation of 37 ± 10 HU could be observed for the six 4DCTs and 6 ± 3 HU for the six 4DCT(CT)s. This resulted in mean lung density differences of 40 ± 10 HU between the 4DCT and the 4DCT(CT) in full inhalation, with differences up to 55 HU being observed.

An overlay of a DRR generated for the EI state of the 4DCT and 4DCT(CT) for CT1–CT6 is shown on the left in Figure 5a–f. The appearance of both DRRs is very similar for all studied patient geometries. Some differences are visible at the boundary of the ribcage, especially for CT4 (Figure 5d) and CT6 (Figure 5f) due to the moving ribcage in the 4DCT, as well as some differences in the rib positions. The relative change of the DRR image intensities from the 4DCT(CT) relative to the values from the DRR from the 4DCT is depicted on the right in Figure 5a–f. Again, some differences are visible due to the moving ribcage and body surface of the 4DCT. As previously explained, the 4DCT(CT) tends to overestimate the lung density for other states than the EE reference state. This is also represented in the relative DRR differences by the blue shading of the lung area. The relative change inside the body, however, is below 5% for >81% of the pixels for all six cases.

FIGURE 5

Overlay of end inhalation digitally reconstructed radiograph (DRR) from 4DCT (green) and 4DCT(CT) (pink) for CT1–CT6 (left) and relative change of the pixel values of the 4DCT(CT) DRR relative to the 4DCT DRR (right)

We compared the original 4DCT to the synthetic 4DCT(CT) and the p4DCT for CT1–CT6. The top row of Figure 6 shows an EI slice of the 4DCT (green) overlaid with the 4DCT(CT) (pink) for all six cases. Due to the static ribcage of the 4DCT(CT), differences are visible in the body surface, ribs, and the anterior edges of the lung. However, for all studied cases, the AP lung expansion is very small and only around 3 mm for CT6, which shows the largest ribcage motion. Due to the DIR uncertainty and interpolation, there are also some differences inside the lung. However, the differences in the anatomy with considerable respiratory changes, such as the diaphragm, are only marginal. For the abdominal organs, the DVF obtained by extrapolation of the DVF of the lungs can lead to some differences (see, e.g. air in Figure 6b). The bottom of Figure 6 shows an overlay of the 4DCT (green) and the p4DCT (pink). The same differences are visible in the ribs, body surface, and at the anterior edge of the lung again due to the static ribcage of the p4DCT. There are no differences in the lung and abdominal organs. An animation of the original 4DCT, 4DCT(CT), and p4DCT of CT2 is given in Figure S4.

FIGURE 6

Overlay of sagittal slice of 4DCT (green) and 4DCT(CT) (pink) in (a–f) and patched 4DCT (pink) in (g–l) for CT1–CT6. All slices show the end inhalation phase and the tumor is visible in all slices

The static dose distribution and the 4D dose distributions calculated on the 4DCT, 4DCT(CT), and p4DCT are shown for the example CT6 in Figure 7a–d. Figure 7e shows the 4D dose differences of the original 4DCT and the 4DCT(CT) for CT6. The dose difference of the p4DCT and the 4DCT(CT) in Figure 7f shows the impact of the density differences in the lung and the DVF interpolation. There is more dose distal to the tumor for the p4DCT due to the higher density in the lung in the 4DCT(CT). On the other hand, Figure 7g shows the impact of the ribcage motion in the dose difference of the original 4DCT and the p4DCT. The results of a 2%/2 mm gamma analysis between the six 4DCTs, the p4DCTs, and the 4DCT(CT)s are listed in Table 1. The gamma pass rates between the 4DCT and 4DCT(CT) are only slightly smaller than the gamma pass rates between the 4DCT and p4DCT. CT3, with the largest tumor isocenter motion (2.8 cm in SI direction), shows larger differences between the three different 4DCT data sets. Another explanation for this, are the larger tumor deformations, which can be seen in Figure 6. V 95% in the CTV and the mean dose in the lung half with the tumor are also listed in Table 1. A clear degradation of the dose coverage in terms of V 95% is visible due to the breathing motion when compared to the static plan, but the dose coverage and the mean lung dose for the 4D plans are very similar among the three 4DCT scenarios, showing that comparable dosimetric indices can be achieved using the 4DCT(CT).

FIGURE 7

Dose distributions and dose differences for CT6 from 4D dose calculation on the original 4DCT (b), 4DCT(CT) (c), and patched 4DCT (d). The static dose distribution is shown in (a) and the chosen field directions are indicated with white arrows. Dose distributions and differences are shown as percentages of the prescribed dose and the same colormaps were used for all dose distributions and dose differences, respectively

TABLE 1

For each of 4DCTs the 2%/2 mm gamma pass rates between the 4DCT and the patched 4DCT (p4DCT), the 4DCT, and the 4DCT(CT) as well as the p4DCT and the 4DCT(CT) are listed

	Gamma pass rate (%)			V 95% CTV (%)				Mean dose healthy lung (%)
	4DCT/p4DCT	4DCT/4DCT(CT)	p4DCT/4DCT(CT)	Static	4DCT	p4DCT	4DCT(CT)	Static	4DCT	p4DCT	4DCT(CT)
CT1	95.0	93.8	98.6	82.8	34.2	33.7	33.5	3.6	3.4	3.5	3.4
CT2	92.9	90.8	97.9	97.5	39.0	38.3	40.1	5.3	5.7	5.7	5.6
CT3	82.5	77.3	94.9	89.9	26.3	27.0	32.3	3.8	4.2	4.2	4.0
CT4	97.8	97.2	99.2	91.9	36.6	37.2	46.4	1.9	1.9	1.9	1.9
CT5	92.6	88.7	96.8	93.5	66.3	67.0	67.8	3.3	3.5	3.5	3.1
CT6	98.9	96.8	97.9	92.6	66.2	66.5	64.7	12.0	12.0	12.0	11.6
Average	93.3	90.8	97.5	91.4	44.7	44.9	47.5	5.0	5.1	5.1	4.9

Note: V 95% and the mean dose to the surrounding healthy lung half were calculated for the static plan and the 4D dose calculations for the 4DCT, p4DCT, and 4DCT(CT).

Abbreviation: CTV, clinical target volume.

DISCUSSION

As shown in this paper, our numerical 4D phantoms are based on a large variety of realistic patient geometries (CTs) and a large motion library based on volunteer 4DMRIs, thus covering a wide variety of 4D treatment scenarios. The phantoms provide an important foundation for further progress in 4D treatment planning and optimization, for the comparison of different motion mitigation techniques under different scenarios as well as for the simulation and comparison of time‐resolved image guidance methods for both proton and photon radiation therapies. The DVFs have been used directly to warp the reference CTs and thus they represent the true underlying motion of each 4DCT(MRI). Due to the availability of these “ground truth” DVFs, the phantoms can also be used for the development or comparison of different DIR algorithms. The use of motion from a 4DMRI instead of a 4DCT allows studying irregular breathing patterns. As 4DMRI leads to no additional dose, long acquisitions are possible, which can for example show the drift of the abdominal organs while the patient lies on the treatment couch. The selected breathing cycles allow for the study of varying degrees of motion irregularity, which is usually not possible with clinically used single‐cycle 4DCTs. Based on this 4D phantom development, we plan several follow‐up studies investigating the impact of motion irregularity on proton treatments. We have also included some proton dose distributions under irregular motion scenarios in Figure S5. To allow others to animate their CTs with different motion patterns, we have made the code for the 4DCT(MRI) generation, together with the complete 4DMRI motion library openly available.

In the original approach developed at our institute for the liver, correspondence between cases was established by generating a mesh from manual landmarks in predefined anatomical locations, leading to a lot of manual work. The approach presented in this paper for lung uses DIR of binary masks to establish correspondence and thus gets rid of this labor intensive step. The workflow could further be automated by automatic segmentation of the lung. A fast, mostly automatic workflow is highly advantageous for adding further CTs and motion patterns to further improve the variability of the phantoms. Additionally, the approach presented here could also more easily be extended to further organs, like the liver, heart or pancreas, whereas motion modeling could be used to further expand the motion library. , , ,

In the literature, a DIR‐based approach to map respiratory motion from 4DMRI onto CTs for abdominal tumors treated with carbon ion therapy has been presented, in which each respiratory phase of the single‐breathing‐cycle 4DMRI is registered to a reference phase of the single‐cycle 4DMRI. They then used multimodality DIR to register a CT and the reference state of the 4DMRI acquired for the same patient to establish correspondence. This approach also assumes a rigid ribcage. Our approach, as described in this paper, together with our previous approach for liver allows the use of a CT and 4DMRI of a different patient or volunteer to be used, which makes it possible to animate one CT with different breathing patterns and consequently create a more extensive data set. Also, our experience shows that the multimodality registration of a CT and an MRI of a different subject is rather challenging due to the different anatomical features. We have therefore applied DIR both to extract the motion from the 4DMRI and to register the binary lung masks to establish correspondence. It is, however, well known that different DIR algorithms can lead to considerably different deformation fields for both intra‐ and interfraction changes imaged for the same patient. , In addition, the registration of different breathing phases from a 4DMRI or a 4DCT is challenged by the presence of sliding organ boundary motion which can lead to discontinuous DVFs. Smooth registration algorithms (like B‐spline) can thus lead to incorrect displacement fields near the boundaries of sliding organs. We use static masks leading to a static ribcage in our phantoms to overcome these potential errors. Further, the registration approach used for the 4DMRIs used a sliding organ mask to deal with this issue. There are also several dedicated registration approaches dealing with sliding boundary motion. ,

A limitation of the phantoms are the density changes in the lung occurring during breathing, which are not considered in the DIR and image warping of the reference CT. Using an EE CT leads to an overestimation of the lung density in the other phases of the synthetic CT. In contrast, the lung density is underestimated for a DIBH CT. The density fluctuations in the lung could be considered by scaling the overall lung density according to observed patient data or by scaling the individual voxel intensities according to the expected density fluctuations calculated from the Jacobian determinant of the DVF. , , , It will be considered in an upgrade of this phantom.

Another limitation is the static ribcage of our phantom. Even for patients breathing shallowly or wearing an abdominal mask restricting the breathing motion, the ribcage is never completely static, and there is also a substantial number of patients undergoing radiotherapy that have non‐shallow breathing patterns. However, for the studied cases using 4DCT motion, we saw a larger impact due to the density differences in the lung than just due to the static ribcage. Additionally, although the provided phantoms may not reflect the full truth, they still provide a number of clinically plausible and realistic scenarios, which is the goal of every numerical phantom. We plan to include both a way to consider the density fluctuations as well as a moving ribcage in the future. As, however, the 4DMRI data used in this project were reconstructed with a method assuming a static ribcage, these data are not suitable for extracting ribcage motion and such motion would have to be gained from 4DCTs or another 4DMRI approach.

We also have to remark on the fact that the motion used for the phantoms was extracted from 4DMRIs of healthy volunteers. While the breathing motion of a lung cancer patient could show different characteristics and especially the motion of cancerous lung tissue could differ from healthy lung tissue leading to differences in tumor motion and deformation, the workflow described in this paper could just as easily be used for 4DMRI acquired for patients for which a CT/4DCT was also acquired. Additionally, subject‐specific 4DCT(MRI)s could also be created using a “motion of the day,” extracted either from a daily 4DMRI or reconstructed from surrogate data, such as from ultrasound.

Whereas in the original XCAT phantom, each organ was modeled as being homogeneous, a lot of effort has been put into a realistic representation of the heterogeneities in the lung by modeling the structures in the lung. , On the other hand, our 4D numerical lung phantoms are completely based on real patient image data and thus directly provide a realistic representation of tissue heterogeneity, which is essential for accurate dose calculation and realistic image guidance simulation, especially for proton therapy.

CONCLUSION

We have developed and validated a method for the generation of 4DCT(MRI) numerical lung phantoms, which provide a realistic representation of a large selection of anatomical and motion scenarios, as well as the accessible deformation vector fields of each 4DCT(MRI). The flexible workflow allows including additional image data and extension of the phantoms to further organs, such as the liver or pancreas. Such phantoms have a very wide range of applications for the investigation of motion mitigation techniques for proton therapy as well as the simulation and comparison of different image guidance techniques.

CONFLICT OF INTEREST

The authors have no conflicts to disclose.

Supporting Information

Click here for additional data file.