If you don't remember your password, you can reset it by entering your email address and clicking the Reset Password button. You will then receive an email that contains a secure link for resetting your password
If the address matches a valid account an email will be sent to __email__ with instructions for resetting your password
GSI Helmholtzzentrum für Schwerionenforschung GmbH, Darmstadt, GermanyInstitute of Condensed Matter Physics, Technical University of Darmstadt, Germany
GSI Helmholtzzentrum für Schwerionenforschung GmbH, Darmstadt, GermanyDepartment of Physics and Astronomy, Louisiana State University, Baton Rouge, LA, USA
Implemented RBE-weighted 4D particle dose alculation for irregular motion.
•
Enables dose calculations on arbitrarily long series of CT images.
•
Impact of irregular motion studied on simulated human phantom CT sequences.
•
Experimental validation of dose calculation algorithm with moving IC array detector.
•
Will enable exploration of effects and study of motion mitigation strategies.
Abstract
Purpose: Highly conformal scanned Carbon Ion Radiotherapy (CIRT) might permit dose escalation and improved local control in advanced stage thoracic tumors, but is challenged by target motion. Dose calculation algorithms typically assume a periodically repeating, regular motion. To assess the effect of realistic, irregular motion, new algorithms of validated accuracy are needed.
Methods:We extended an in-house treatment planning system to calculate RBE-weighted dose distributions in CIRT on non-periodic CT image sequences. Dosimetric accuracy was validated experimentally on a moving, time-resolved ionization chamber array. Log-file based dose reconstructions were compared by gamma analysis and correlation to measurements at every intermediate detector frame during delivery. The impact of irregular motion on treatment quality was simulated on a virtual 4DCT thorax phantom. Periodic motion was compared to motion with varying amplitude and period baseline drift. Rescanning as a mitigation strategy was assessed on all scenarios.
Results:In experimental validation, average gamma pass rates were 99.89+-0.30% for 3%/3 mm and 88.2+-2.2% for 2%/2 mm criteria. Average correlation for integral dose distributions was 0.9900.002. Median correlation for single 200 ms frames was 0.9470.006. In the simulations, irregular motion deteriorated V95 target coverage to 81.2%, 76.6% and 79.0% for regular, irregular motion and irregular motion with base-line drift, respectively. Rescanning restored V95 to 98% for both scenarios without baseline drift, but not with additional baseline drift at 83.7%.
Conclusions: The validated algorithm permits to study the effects of irregular motion and to develop and adapt appropriate motion mitigation techniques.
Intrafractional baseline shift or drift of lung tumor motion during gated radiation therapy with a real-time tumor-tracking system a preliminary version of this study was presented at the 55th Annual Meeting of the American Society for Radiation Oncology.
]. Therefore, treatment methods have been developed in attempt to handle tumor motion during radiotherapy including gating, ITVs, rescanning, tumor tracking [
Pfeiler T, Ahmad Khalil D, Bäumer C, Blanck O, Chan M, Engwall E, et al. 4D robust optimization in pencil beam scanning proton therapy for hepatocellular carcinoma. In J Phys Conf Ser, volume 1154. Institute of Physics Publishing. ISSN 17426596, 012021. doi: 10.1088/1742-6596/1154/1/012021.
Delivery quality assessment is an integral part of clinical radiotherapy. Typically, treatment planning, simulations and dose reconstructions are performed to simulate expected dose distributions and assess delivered dose distributions respectively [
Log file-based dose reconstruction and accumulation for 4D adaptive pencil beam scanned proton therapy in a clinical treatment planning system: Implementation and proof-of-concept.
]. However, only periodic repetition of the respiratory motion is currently assumed for treatment planning on a time-resolved computed tomography (4DCT) [
]. 4D dose calculations take into account the inter- and intrafractional motion variability of recorded motion signals when assigning motion phases, but still assumes that one 4DCT would be representative for an entire treatment fraction [
] reported interfractional variations in the motion amplitude of more than 5 mm in nearly half of the patients with average tumor motion amplitudes greater than 7 mm and intrafractional changes in the motion amplitudes of up to 10 mm, increasing with total motion amplitude. Typical baseline drifts of greater than 5 mm in cranio-caudal direction were reported over the course of treatment for lung tumors and baseline drifts greater than 2 mm were reported for liver lesions. Consequently, dose calculation algorithms, which are able to handle complex motion patterns, are needed to investigate and accurately assess treatments with various motion mitigation techniques and to determine adequate safety margins during planning.
Current online imaging and motion monitoring methods are not sufficient to verify 4D dose calculations and do not measure directly stopping powers along the entire beam path. 4DCT images cannot be acquired for the time span of an entire treatment delivery, due to the unacceptably high doses and required imaging frequency [
Integrated radiotherapy imaging system (IRIS): design considerations of tumour tracking with linac gantry-mounted diagnostic x-ray systems with flat-panel detectors.
]. Online images, such as cine MRI and cone beam CT images, do typically not include deformations of the surrounding tissues directly. In order to deduce the resulting changes of the penetration depth of the ion beam, the obtained images have to be translated into stopping powers by registering them to a CT image. To overcome this problem, Phillips et al. [
] proposed a method that deduces water equivalent path lengths (WEPLs) from a single CT and translates them into a beam-specific WEPL space. Tumor motion is then simplified as a translation in the WEPL space and deduced from surrogate data. Deformations of the surrounding tissue are neglected. Meschini et al. applied this method to RBE-weighted dose calculations for carbon beams [
]. The reduction of needed 4D imaging of the tumor makes this method applicable only for dose calculations based on 2D information about the tumor motion. That data can be recorded during treatment, which makes this method clinically applicable. However, this method is limited in case of strong soft tissue deformations.
Full information about geometry changes is necessary and already available in two kinds of application: First, simulation studies can be used to identify tumor locations with critical dose degradation due to irregular motion. The required long image sequences can be created by using deformable image registration based on one reference 3DCT and time-resolved magnetic resonance imaging (4D-MRI) data [
], which was used in this study. Second, testing and benchmarking of new motion mitigation techniques can be done in simulations and experiments with simplified and anthropomorphic phantoms [
Ehrbar S, Perrin R, Peroni M, Bernatowicz K, Parkel T, Pytko I, et al. Respiratory motion-management in stereotactic body radiation therapy for lung cancer – A dosimetric comparison in an anthropomorphic lung phantom (LuCa). Radiother Oncol 2016;121. doi:10.1016/j.radonc.2016.10.011.
]. In this case, detailed information about the phantom motion is known and virtual CTs can be generated.
For densly ionizing irradiation, the biological effect on irradiated cells is higher than for a photon irradiation with the same dose. The concept of relative biological effectiveness (RBE) and RBE-weighted doses links absorbed ion and photon doses with the same biological endpoint. For protons, a constant cell killing RBE of 1.1 can be assumed [
Relative biological effectiveness (RBE) values for proton beam therapy. Variations as a function of biological endpoint, dose, and linear energy transfer.
]. For carbon beams, in contrast, the RBE depends on many factors as the particle spectrum and local energy transfer (LET) and is non-linear with dose [
Quantification of the relative biological effectiveness for ion beam radiotherapy: Direct experimental comparison of proton and carbon ion beams and a novel approach for treatment planning.
]. Consequently, dose warping from on phase CT to a reference CT cannot be applied. Instead, the contributions of single pencil beams need to be added properly in the reference CT.
The objective of this study was to provide an experimentally validated dose calculation algorithm that can calculate absorbed and RBE-weighted dose distributions on arbitrarily long sequences of CT images, as well as giving an example of the workflow to apply this algorithm in simulation studies.
2. Materials and methods
In order to be able to calculate doses on an arbitrarily long sequence of CT images, the treatment planning system TRiP98 developed at GSI [
Zeitlich aufgelöste Dosisberechnung für gescannte Kohlenstoffionen auf sequentiellen Bilddaten zur Erfassung unregelmäßiger Bewegungen (Master Thesis).
] was modified. The changes are described in section 2.1. We experimentally validated the calculations of absorbed dose distributions for irregular motion, described in section 2.2,3.1. Finally, we applied the algorithm for RBE-weighted dose calculations to a virtual phantom as described in 2.3, 3.2.
] was designed for the calculation of biologically effective doses in particle therapy on 4DCT imaging data, which were assumed to be periodically repeating. The required data therefore consist of the 4DCT and corresponding transformation maps relating the several motion states to one reference state. This limited amount of data can be read at once and held in memory during the entire dose calculation. In the case of irregular motion, CT sequences covering the entire course of treatment can be substantially larger. In order to handle the additional amount of data, the 4D-dose calculation algorithm was altered to sequentially calculate the dose contributions of the motion states and only hold the data necessary for the current motion state in memory.
In this section, the order of calculation steps in the existing dose calculating algorithm is outlined. Then, the changes made to the algorithm are detailed.
The relative biological effectiveness is defined as the ratio of a photon dose and a particle dose yielding the same biological effect. In order to calculate the biological effect of an irradiation with ion irradiation, the dose dependent cell survival S is estimated with a linear-quadratic model (LQM) [
For the 4D dose calculation, an irradiation of points on a 4-dimensional grid is assumed. The first coordinates of a spot, labeled by index i, is the motion state s during which it was irradiated. The other coordinates are the position and in the iso-energetic slice (IES), which corresponds to the initial beam energy . In order to take into account the entire particle spectrum, the dose averaged coefficients
(2)
(3)
with the normalization factor
(4)
are calculated based on the local effect model (LEM) [
], the virtual index l of the beam component comprises the irradiated raster spot i, particle type T, energy E and the water equivalent depth of the beginning and end of the corresponding voxel in the body. The weight represents the number of particles with these properties that contributes to the dose. The particle spectra in T and E and the resulting contributions and were precomputed for available initial beam energies and water equivalent depths k and stored in a data base. This reduces the calculation of the dose averaged coefficients to
(5)
(6)
(7)
Here, is the weight in a linear interpolation for the voxel between precomputed water equivalent depths before and behind the actual depth. In the standard algorithm for periodic motions, computation proceeds voxel by voxel on one CT defined as reference. The loop over the motion states is nested inside the outer voxel loop. The changes in voxel position between the reference CT and other motion state CTs are looked up from a precalculated transformation map. The sums in Eq. (5), (6), are updated by calculating their addends on the corresponding state CTs. Finally, the normalization of Eq. 5, is performed. The physically absorbed dose in the reference state CT voxel is given by
The RBE-weighted dose is then determined by finding the X-ray dose with the same survival S as the one calculated for ions by inverting Eq. 1:
(10)
Final results are stored into a grid with the same dimension as the CT image and then, computation proceeds for the next reference CT voxel.
To perform dose calculations considering irregular motion, long series of CTs and transformation maps need to be handled. This is achieved by treating the motion states sequentially. This means that the loop over the motion states is changed to be the outer-most one. Inside each motion state, all possible calculations proceed voxel by voxel. In the following, the necessary changes to the algorithm are described. The normalization of Eq. 5, can only be performed when the treatment was completed and, thus, the sum in Eq. 7 was completed too. In order to collect the contributions to those three equations, additional three CT sized grids are introduced. To reduce the amount of required working memory for the motion state loop, only the CT and transformation map for the currently considered motion state are held in storage. The contributions to Eq. (5), (6) (without normalization) were calculated for all points irradiated during this motion state, mapped to the corresponding position in the reference state and added to the corresponding grids. After iterating over all motion states, the normalization of Eq. 5 is performed and physical and RBE-weighted dose distributions in the reference state are computed using equations Eq. (8), (9). In order to study the dynamics of an entire delivery, an output of results for single motion states was implemented. For this, one slice of the reference CT can be selected. The dose absorbed in this slice during each single motion state is stored. These slices were then appended to form a new CT where the slice index indicates the motion state. In the same manner, an integral slice, showing the physical dose that was accumulated up to the indexed motion state, was also stored in order to reconstruct the build-up of the absorbed dose over time.
2.2 Experimental validation
For the experimental validation, the goal was to test if the modified dose calculation algorithm is able to reproduce inhomogeneous dose distributions caused by the interplay between scanned beams and irregular motion patterns. Experiments were performed at the horizontal beam line in treatment room 1 at CNAO, Pavia, Italy. The experimental setup is schematically displayed in Fig. 1.A 2D ionization chamber (IC) array detector (PTW Octavius 1500XDR; PTW, Freiburg, Germany) was mounted on a programmable linear stage (PI M-414.2PD, controller PI C-884.4DC; Physik Instrumente, Karlsruhe, Germany). The detector consists of 4.4 4.4 x 3mm3 big IC chambers in a chessboard like arrangement. Measurements were performed in a time resolved mode with a frame length of 200 ms. It is placed inside a dedicated holder behind 5 mm of PMMA. The linear stage was programmed to move in sinusoidal patterns of the form
(11)
with amplitude and period variations between subsequent periods, each period ranging from one minimum to another. Here, and denote the motion amplitude and period of the current motion and the time of the precedent minimum. The average period T was set to 5 s. Motion periods and amplitudes were varied to uniformly distributed random values in the range ] and , respectively. A series of seven measurements labeled I to VII with different motion patterns was performed. The motion settings are summarized in Table 1. 238.63 Me V/u scanned carbon ion beams with particles per spot were delivered across a 60 60 mm2 grid. Spot spacing was 2 mm and the beam was assumed to be Gaussian shaped with a full width half maximum of 5.78 mm. The beam was scanning line by line with alternating direction over the target volume. The motion of the detector was measured with an optical distance laser sensor (SICK OD100-35P840; SICK AG, Waldkirch, Germany).
Each detector time frame was defined as one motion state. The position and number of particles for the beam spots irradiated during one time frame were retrieved from the delivery system log files [
] by a LabVIEW based software developed in house. This delivery information was extracted from the magnet currents of the scanning magnets and the signal from the IC intensity monitor integrated in the nozzle. To synchronize the clocks of the IC array detector and the dose delivery system (DDS), the correlation coefficient between the signal recorded by the intensity monitor and the recorded energy deposition in the IC array detector during the same detector time frames was maximized. To do so, a first estimate of the clock offset was deduced by matching the end of the last frame with measured dose and the timestamp of the last delivered beam spot. Then the correlation was sampled by scanning possible values for the clock offset in range of 2000 ms centered around this estimate with a resolution of 100 ms. The hardware counts of the intensity monitor were logged by the DDS at every progression to the next raster spot or at least every 50 ms. At the transition between two subsequent detector time frames, the IC counts of the corresponding log entry were split between those time frames according to the time fraction in each frame. In a second step, the scan range was reduced to the range between the sampled value before and after the maximum and the correlation was sampled again, but with a ten fold smaller step size. This procedure was repeated once more to get a final step size of 1 ms. The average detector displacement during one detector frame was used as the translation vector for dose calculations.
The log-files were converted to TRiP98 file format. The absorbed dose distribution was calculated using the new algorithm inside a 1 mm grid in water. The physical dose distribution was calculated and compared to the measurement of the IC array detector in a water equivalent depth of 13.3 mm. This is the sum of the water equivalent thicknesses of the PMMA holder (WEPL = 1.15 [
] with local dose criterion was used to quantify the agreement of measurement and dose reconstruction for the expected inhomogeneous dose distribution. The measured dose in each IC was used as reference dose to determine if this IC passes the gamma test. Acceptance criteria of 3%/3 mm and of 2%/2 mm were selected. Below a dose threshold of , corresponding to 5% of the planned dose of 3Gy that would have been delivered to the homogeneous target region of a static delivery, an absolute dose acceptance criterion of was used. In addition, Pearson’s correlation coefficient R between the dose values in the dose reconstruction and the measured IC dose values was calculated. For both comparisons, only IC detectors with a non-zero dose were considered. The evolution of R over the course of the treatment was calculated for the dose distributions inside single detector time frames (in the following called differential dose distributions) as well as with intermediate integral doses up to a given motion state. In order to determine the number of motion states per period required for precise dose reconstructions, the entire dose reconstruction procedure was repeated, binning subsequent detector frames in various multiples to one motion state and the pass rate of the final dose distribution was compared dependent on the time resolution.
2.3 Dose calculation for irregular motion using a virtual phantom
In this section, a demonstration on how the dose calculation algorithm can be used in combination with a virtual phantom to study the effects of irregular motion on the quality of patient treatments are presented.
As a preliminary test of the ability of the dose calculation algorithm, a limited set of virtual CTs was created using the software XCAT [
]. From this set, a long sequence of images depicting non-periodic motions was created. The procedure is described in the following:
A spherical tumor with a diameter of 40 mm was placed inside the right lung of a phantom (see Fig. 2). 25 virtual CTs depicting half a respiratory cycle from end-exhale (state 0) to end-inhale (state 24) were created. As parameters, a motion amplitude of the diaphragm’s dome of 40 mm and an anterior-posterior motion of the rib cage with an amplitude of 20 mm were set and both displacements were set to increase linearly with the state index. The used diaphragm and anterior-posterior motion amplitudes are twice as big as the values recommended for an average breathing cycle. For the respiratory motion during treatment, a motion pattern as proposed by Lujan et al. [
] was selected. During breathing cycle k, so for , the displacement d of the diaphragm’s dome from its standard end-exhale position in cranial direction was given by the equation
(12)
Fig. 2Virtual CT containing a spherical tumor created using XCAT. The image shows a superposition of the motion states 7 and 19, which correspond the end-exhale and inhale states of the planned motion. The arrows are used to visualize the corresponding transformation map.
Here, and represent motion amplitude of the diaphragm’s dome and period for respiratory cycle k, respectively. is the start position of the first breathing cycle and is the velocity of a linear baseline drift. Three motion scenarios were considered for the dose calculation:
Periodic: Periodic motion from Eq. 12 with and for all periods and no baseline drifts.
Fixed baseline: Nonperiodic motion, with and randomly generated with a uniform distribution between +-25% around and , and no baseline drift.
Baseline drift: The same motion as in the fixed baseline scenario, but with an additional linear baseline drift of .
Treatment delivery progression was simulated for the GSI facility with a constant beam intensity during a spill for single scan deliveries and 20 slice-by-slice rescan deliveries. For the deliveries with a single scan, the irradiation lasted 124 s for field 1 and 128 s for field 2. For deliveries with 20 rescans, the irradiations lasted 804 s for field 1 and 642 s for field 2. Each breathing period was divided into ten motion states of equal length. For each state, the center-of-mass of the displacement of the diaphragm was computed. An offset was chosen such that for the baseline drift motion scenario and 20 rescans, the smallest diaphragm displacement corresponds to the simulated end-exhale CT with index 0. All motion patterns started at , corresponding to CT number 7. For the dose calculation, this motion state was defined the reference state . The mean diaphragm displacement for each motion state j was calculated and the precomputed state with the closest diaphragm displacement was selected. A series of more than 1000 CTs was created that is a series of selections of the 25 virtual CTs created before.
In addition to the CT images, the corresponding transformation maps, that link the state CTs to the reference CT, are needed. The forward transformation vector fields for were directly exported from XCAT. For each voxel in the reference state , they define a vector to the corresponding position in the other states i. For , the transformation was obtained by mapping from state to the repetition of state i in a second, identical breathing cycle, . In order to get the inverted transformation fields, XCAT was rerun for each motion state i. The state i was assigned to be the start phase and the first repetition of the reference state was the end state. Now the vector field was exported. The output data produced by XCAT are the attenuation in binary format for the virtual CT images and vector fields in a text format for the transformation maps. Those files were converted into the corresponding input formats for TRiP98 using a software developed in house. The series of transformation files for the dose calculation was created with the same map as for the corresponding CT images.
For the purpose of treatment planning, the motion was assumed to be periodic. The motion used during treatment planning was split into ten phases of equal length with a mean displacement corresponding to the simulated CTs and 7. This corresponds to a maximum tumor displacement of 9.83 mm in cranio-caudal direction, −3.82 mm in antero-posterior direction and −0.18 mm in left–right direction. The union of the tumor contours in these 10 states was defined as the ITV. Two opposed fields were optimized simultaneously for a homogeneous RBE-weighted dose of 8 Gy in the ITV. The same motion pattern was assumed for both fields. This high dose was selected due to encouraging results in recent dose escalation studies with carbon ions for non-small cell lung cancer [
]. The resulting dose distributions were compared by calculating the DVH metrics such as and on the gross tumor volume. is the fraction of the volume which receives a dose of the prescribed dose and quantifies target coverage. and denote the dose that is exceeded in 5% and 95% of the target volume, respectively. The difference quantifies the width of the dose fall-off in a dose volume histogram around the prescription dose. The smaller it is, the more homogeneous is the dose distribution with a mean dose close to prescription.
3. Results
3.1 Experimental validation
We provide detailed results for each analysis step of measurement VII with motion amplitude and period both varying by 50% as this is the most complex motion. Final results for all other measurements are listed in Table 2.
Table 2Comparison of dose reconstruction and measured dose distribution for measurements with different average motion amplitude A and variabilities of motion amplitude and period T. The average motion period was always set to 5 s.
The recorded motion for VII is displayed in Fig. 3. Fig. 4 shows the correlation between the measured counts of the beam monitor and the integral dose in the single ICs for each matched frame of the IC detector array. A squared correlation coefficient of was obtained. The reconstructed dose is compared with the measurement in Fig. 5. The simulated dose was convoluted with the size of the single IC detectors in the array. In Fig. 5(a) an inset of reconstructed and measured dose is shown. In Fig. 5(b) the corresponding correlation plot of measured vs. reconstructed dose is shown. A linear regression yields a correlation coefficient . The correlation coefficient was also determined for single time frames of the ionization chamber array detector. In Fig. 6, the evolution of Pearson’s R over the course of the delivery is displayed. The green dots show the correlation coefficient of the dose only in the single time frames of the detector. The black line shows how the correlation coefficient of the integral dose and the simulation evolves. The results for the integral dose distributions of all measurements are summarized in Table 2.
Fig. 3Recorded motion of the ionization chamber array detector for measurement VII. Black dots represent measured positions, the red line is an interpolation.
Fig. 4Scatter plot showing the measured number of counts in the intensity monitor over the detector response of the IC array detector for each 200 ms time frame for measurement VII. The red line shows an unconstrained linear fit and was added to guide the eye.
Fig. 5Comparison of reconstructed and measured dose distribution for measurement VII. (a) The small homogeneous squares in the foreground show the measured dose values in the single ionization chambers (IC). The IC contour size was reduced for better visibility. In the background the reconstructed dose distribution is displayed with 1 mm resolution. (b) Correlation of measured and reconstructed dose distribution for all ICs with measured dose. The red linear fit was added to guide the eye.
Fig. 6Left: Time resolved evolution of the correlation coefficient for measurement VII. Right: Binned projection of the correlation coefficients of the differential measurements on the y axis.
The evolution of the correlation coefficient between measured and simulated integral dose distribution is displayed in Fig. 7 for all seven cases. The correlation is stable above 0.9 and converges quickly. All measurements show a common evolution. In Fig. 8, the distributions of correlation coefficients for the differential dose distributions in the single detector frames is shown. Detector frames without dose, which occurred for example during spill pauses, were excluded from the selection. The median of the correlation coefficient for each measurement is listed in Table 2. The dependence of the pass rate on the length of each motion state is depicted in Fig. 9. It can be seen that for a average motion amplitude of 20 mm, a frame duration of 400 ms still yielded a pass rate greater than 99% in all cases. In motion scenario III with an amplitude of 30 mm, in contrast, the pass rate drops already for this frame rate to 97.84%. For motion I with an amplitude of 10 mm, the pass rate remained 100% up to a frame duration of 800 ms.
Fig. 7Evolution of the correlation between the measured and reconstructed integral dose distributions over the course of the measurements listed in Table 2.
Fig. 8Distribution of correlation coefficients between reconstructed and measured differential dose distribution for single detector frames for all seven measurements. The histograms were created the same way as the one on the right of Fig. 6 and are shown as violin plots overlaid with box plots. The horizontal lines depict the median and the boxes the interquartile range of t.he distribution.
3.2 Dose calculation for irregular motion using a virtual phantom
The simulated distributions of RBE-weighted dose for the motion scenarios described in section 2.3 are displayed in Fig. 10. The resulting DVH metrics for and are listed in Table 3. It can be seen, that for the short single deliveries, the interplay effect causes inhomogeneities of similar size. Rescanning improves dose homogeneity and target coverage in all three scenarios, but for the linear baseline drift scenario, the increased treatment duration causes the tumor to move significantly out of the target volume, causing a lower target coverage of (see Fig. 10f).
Fig. 10Calculated RBE-weighted dose distributions for the three motion scenarios described in section 2.3 (a), (c), and (e) for a single scan and (b), (d), and (f) for 20 rescans. The black circle marks the contour of the spherical lesion.
In this work, we presented a dose calculation algorithm to calculate RBE-weighted doses for complex motion patterns on arbitrarily long sequences of CT images. This algorithm was validated experimentally by performing logfile-based reconstructions of the absorbed dose distributions for seven different motion scenarios and comparing them to measurements with an IC array detector. Excellent agreement between dose reconstructions and measurements was achieved for integral doses as well as for the dynamics of the dose delivery down to the level of single 200 ms time frames. It was demonstrated that the algorithm can be applied for the calculation of RBE-weighted doses by using images from a virtual phantom.
The experimental validation of the dose calculation showed, that the modifications necessary for dose calculations with complex motion patterns cause no relevant additional uncertainties compared to the ones already present for periodic motion. The experimental data strongly supports the reliability of simulation studies performed with the described dose calculation algorithm. In addition, the results demonstrate the high reliability of the methods used to create the log-file based dose reconstructions. In experimental studies, this allows for examination of the quality of a delivery not only inside a detector plane, but also in volumes, for example calculated on a patient 4DCT. Nevertheless, it is recommended to use a detector to verify the quality of the dose reconstruction in at least one plane if possible. This is done for example by Lis et al. [
] for the standard algorithm for periodic motion. Another important finding is that for an average motion amplitude of 20 mm, about ten motion states per period are required to guarantee a pass rate greater than 99%, which corresponds to a state duration of between 400 ms and 600 ms in the measurements presented here. As expected, the required time resolution must be finer for bigger motion amplitudes and can be coarser in the case of smaller amplitudes. This result is in perfect agreement with the one presented by Bert et al. [
] and extends it from regular to irregular motion. Deviations between measured and calculated doses were observed in regions with high dose gradients. This is apparent from the high gamma pass rates. Dose differences are caused by the residual motion during calculated motion states that were not accounted for by the dose calculation. In a clinical case, motion mitigation techniques would be necessary for motion amplitudes as big as the ones used during this study. In a motion compensated treatment, dose gradients are small inside the target regions. Therefore, the precision of the algorithm should be sufficient.
The simulations of the RBE-weighted doses demonstrated a workflow how the virtual phantom XCAT can be used for simulation studies in order to test motion mitigation techniques and to identify tumor locations were the dose distribution is sensitive to irregular respiratory motion. In the given example, the technique of rescanned ITVs was able to deliver a homogeneous dose to a tumor even in the case of variable motion amplitude and period with stable baseline. In the case of a linear baseline drift, in contrast, significant underdosage was observed at the edges of the target volume. This finding underlines the need of sufficient safety margins, an active motion mitigation technique or the restriction of the treatment to a stable gating window [
The results of the simulations on the virtual phantom presented in this paper are strongly supported by the experimental validation presented in section 2.2 and section 3.1. The experimental validation was only possible for absorbed doses and was performed only in one water equivalent depth in the plateau region of the Bragg peak. Nevertheless, it incorporates all motion related features that also apply to the calculation of RBE-weighted doses. The TRiP98 dose calculation algorithm was validated experimentally for periodic motions using radiochromic films [
]. As only the management of motion in the dose calculation was altered and the functionality of this component was shown to work in 2D, we assume that the results can be translated to 3D dose distributions and RBE-weighted doses. Nevertheless, more experiments are necessary and will be published in the future. The results of this work will enable to evaluate the quality of dose distributions during the future development of motion mitigation techniques, capable of handling irregular motion. This is true for both, simulations on virtual or real patient imaging data, and in experiments where it gives information about regions not covered by detectors, dose distributions inside a patient instead of in a phantom and about the time structure of the dose delivery. It was shown that the dose calculation algorithm is able to reproduce the final dose distributions as well as the dynamics of the dose build up during treatment. This makes it possible to identify even the effect of short interruptions of regular breathing like coughing.
The experimental validation of the dose calculation algorithm is limited by the simplicity of the geometry, which was selected for this proof-of-principle study. More measurements with complex or anthropomorphic phantoms [
Ehrbar S, Perrin R, Peroni M, Bernatowicz K, Parkel T, Pytko I, et al. Respiratory motion-management in stereotactic body radiation therapy for lung cancer – A dosimetric comparison in an anthropomorphic lung phantom (LuCa). Radiother Oncol 2016;121. doi:10.1016/j.radonc.2016.10.011.
] would be desirable and will be carried out in the future. However a more complex geometry also makes the interpretation of the measurement data more challenging.
The calculation of the biological doses described in 2.3, 3.2 are limited to the purpose of demonstrating the abilities of the dose calculation algorithm described in section 2.1. Systematic studies about the effect of irregular motion during carbon ion therapy need to be performed in the future. A study with real patient 4DCT-MRI data would allow for a more realistic study [
], even though the amount of work is not feasible for clinical practise now. With further automation, 3D imaging data might become available. For current research purposes, in contrast, XCAT was shown to be a versatile tool to create virtual CT data sets for exploratory studies as it allows for the isolated examination of the effect of different motion patterns as uncertainties induced by imaging, changes in patient anatomy and deformable image registration can be eliminated. The problems in the creation of CT data sets using XCAT described in [
] are known to the authors. The post-processing tool presented in this reference will be included into the workflow for later studies. From the simulations, it is apparent that rescanned ITVs are a possible approach to achieve homogeneous dose distributions and good target coverage for the cases of variable motion amplitudes. But this comes with the detriment, that increased treatment times cause a significant impact of baseline drifts. Robust 4D optimization [
]. Here, only a linear baseline drift was considered. The used drift velocity for the tumor is in good agreement with the mean trend reported by Takao et al. [
Intrafractional baseline shift or drift of lung tumor motion during gated radiation therapy with a real-time tumor-tracking system a preliminary version of this study was presented at the 55th Annual Meeting of the American Society for Radiation Oncology.
]. It is expected that higher variation in baseline drifts will further exacerbate the dose degradation and will pose a challenge in the development of future motion mitigation techniques.
5. Conclusions
In this work, we presented an algorithm for 4D dose calculation in carbon ion therapy that is able to handle irregular motion. An experimental validation demonstrated its ability to reproduce measurement data for absorbed doses with high precision. We presented a case study, in which the algorithm was used in combination with a virtual anthropomorphic phantom to perform simulation studies. In the future, this procedure can be used to identify tumor locations with high sensitivity to motion irregularities. The presented algorithm will be used to test motion mitigation techniques as real time tracking [
] in the presence of irregular motion experimentally and in simulations.
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgements
We would like to thank Chistoph Schuy and Ulrich Weber for providing the PTW Octavius 1500XDR detector, and the GSI workshop for their support in building the phantoms for our experimental setups.
Intrafractional baseline shift or drift of lung tumor motion during gated radiation therapy with a real-time tumor-tracking system a preliminary version of this study was presented at the 55th Annual Meeting of the American Society for Radiation Oncology.
Pfeiler T, Ahmad Khalil D, Bäumer C, Blanck O, Chan M, Engwall E, et al. 4D robust optimization in pencil beam scanning proton therapy for hepatocellular carcinoma. In J Phys Conf Ser, volume 1154. Institute of Physics Publishing. ISSN 17426596, 012021. doi: 10.1088/1742-6596/1154/1/012021.
Log file-based dose reconstruction and accumulation for 4D adaptive pencil beam scanned proton therapy in a clinical treatment planning system: Implementation and proof-of-concept.
Integrated radiotherapy imaging system (IRIS): design considerations of tumour tracking with linac gantry-mounted diagnostic x-ray systems with flat-panel detectors.
Ehrbar S, Perrin R, Peroni M, Bernatowicz K, Parkel T, Pytko I, et al. Respiratory motion-management in stereotactic body radiation therapy for lung cancer – A dosimetric comparison in an anthropomorphic lung phantom (LuCa). Radiother Oncol 2016;121. doi:10.1016/j.radonc.2016.10.011.
Relative biological effectiveness (RBE) values for proton beam therapy. Variations as a function of biological endpoint, dose, and linear energy transfer.
Quantification of the relative biological effectiveness for ion beam radiotherapy: Direct experimental comparison of proton and carbon ion beams and a novel approach for treatment planning.
Zeitlich aufgelöste Dosisberechnung für gescannte Kohlenstoffionen auf sequentiellen Bilddaten zur Erfassung unregelmäßiger Bewegungen (Master Thesis).
00324-0&id=gr1.jpg){kind=link}
00324-0&id=gr2.jpg){kind=link}
00324-0&id=gr3.jpg){kind=link}
00324-0&id=gr4.jpg){kind=link}
00324-0&id=gr5.jpg){kind=link}
00324-0&id=gr6.jpg){kind=link}
00324-0&id=gr7.jpg){kind=link}
00324-0&id=gr8.jpg){kind=link}
00324-0&id=gr9.jpg){kind=link}
00324-0&id=gr10.jpg){kind=link}