Advertisement

Extension of RBE-weighted 4D particle dose calculation for non-periodic motion

Published:October 26, 2021DOI:https://doi.org/10.1016/j.ejmp.2021.10.009

      Highlights

      • 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.990±0.002. Median correlation for single 200 ms frames was 0.947±0.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.

      Keywords

      1. Introduction

      Scanned charged particle therapy has become an established therapy for static tumors, offering more conformal dose distributions than photon therapy [
      • Durante M.
      • Loeffler J.S.
      Charged particles in radiation oncology.
      ,
      • Grau C.
      • Durante M.
      • Georg D.
      • Langendijk J.A.
      • Weber D.C.
      Particle therapy in Europe.
      ]. This high degree of conformity poses geometric and dosimetric uncertainties when accounting for inter- and intrafractional motion [
      • Riboldi M.
      • Orecchia R.
      • Baroni G.
      Real-time tumour tracking in particle therapy: technological developments and future perspectives.
      ,
      • Bert C.
      • Durante M.
      Motion in radiotherapy: particle therapy.
      ,
      • He P.
      • Mori S.
      Perturbation analysis of 4d dose distribution for scanned carbon-ion beam radiotherapy.
      ]. Tumors in the thorax and abdomen are strongly affected by respiratory motion and its variability [
      • Takao S.
      • Miyamoto N.
      • Matsuura T.
      • Onimaru R.
      • Katoh N.
      • Inoue T.
      • 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.
      ,
      • Dhont J.
      • Vandemeulebroucke J.
      • Burghelea M.
      • Poels K.
      • Depuydt T.
      • Van Den Begin R.
      • et al.
      The long- and short-term variability of breathing induced tumor motion in lung and liver over the course of a radiotherapy treatment.
      ]. Therefore, treatment methods have been developed in attempt to handle tumor motion during radiotherapy including gating, ITVs, rescanning, tumor tracking [
      • Riboldi M.
      • Orecchia R.
      • Baroni G.
      Real-time tumour tracking in particle therapy: technological developments and future perspectives.
      ,
      • Bert C.
      • Herfarth K.
      Management of organ motion in scanned ion beam therapy.
      ] and 4D optimization [
      • Graeff C.
      • Constantinescu A.
      • Lüchtenborg R.
      • Durante M.
      • Bert C.
      Multigating, a 4D optimized beam tracking in scanned ion beam therapy.
      ], as well as robust optimization [

      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.

      ,
      • Wolf M.
      • Anderle K.
      • Durante M.
      • Graeff C.
      Robust treatment planning with 4D intensity modulated carbon ion therapy for multiple targets in stage IV non-small cell lung cancer.
      ,
      • Kanai T.
      • Paz A.
      • Furuichi W.
      • Liu C.S.
      • He P.
      • Mori S.
      Four-dimensional carbon-ion pencil beam treatment planning comparison between robust optimization and range-adapted internal target volume for respiratory-gated liver and lung treatment.
      ].
      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 [
      • Meijers A.
      • Jakobi A.
      • Stützer K.
      • Guterres Marmitt G.
      • Both S.
      • Langendijk J.A.
      • et al.
      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) [
      • Vedam S.S.
      • Keall P.J.
      • Kini V.R.
      • Mostafavi H.
      • Shukla H.P.
      • Mohan R.
      Acquiring a four-dimensional computed tomography dataset using an external respiratory signal.
      ]. 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 [
      • Meijers A.
      • Knopf A.C.
      • Crijns A.P.
      • Ubbels J.F.
      • Niezink A.G.
      • Langendijk J.A.
      • et al.
      Evaluation of interplay and organ motion effects by means of 4D dose reconstruction and accumulation.
      ]. This means that a variable breathing amplitude can accounted for, but varying amplitudes or baseline drifts cannot. Dhont et al. [
      • Dhont J.
      • Vandemeulebroucke J.
      • Burghelea M.
      • Poels K.
      • Depuydt T.
      • Van Den Begin R.
      • et al.
      The long- and short-term variability of breathing induced tumor motion in lung and liver over the course of a radiotherapy treatment.
      ] 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 [
      • Bert C.
      • Durante M.
      Motion in radiotherapy: particle therapy.
      ], even though techniques to reconstruct 3D CT images from 2D cone beam CT images were presented [
      • Li R.
      • Jia X.
      • Lewis J.H.
      • Gu X.
      • Folkerts M.
      • Men C.
      • et al.
      Real-time volumetric image reconstruction and 3D tumor localization based on a single x-ray projection image for lung cancer radiotherapy.
      ]. Tumor motion itself can be monitored directly using 4D cine MRI [
      • Raaymakers B.W.
      • Raaijmakers A.J.
      • Lagendijk J.J.
      Feasibility of MRI guided proton therapy: magnetic field dose effects.
      ], ultrasound [
      • O’Shea T.
      • Bamber J.
      • Fontanarosa D.
      • Van Der Meer S.
      • Verhaegen F.
      • Harris E.
      Review of ultrasound image guidance in external beam radiotherapy part II: Intra-fraction motion management and novel applications.
      ], and fiducial markers [
      • Berbeco R.I.
      • Jiang S.B.
      • Sharp G.C.
      • Chen G.T.
      • Mostafavi H.
      • Shirato H.
      Integrated radiotherapy imaging system (IRIS): design considerations of tumour tracking with linac gantry-mounted diagnostic x-ray systems with flat-panel detectors.
      ] or indirectly by measuring a surrogate signal with known correlation to the tumor position [
      • Dong B.
      • Graves Y.J.
      • Jia X.
      • Jiang S.B.
      Optimal surface marker locations for tumor motion estimation in lung cancer radiotherapy.
      ]. Reviews about imaging techniques in particle therapy can be found in the literature [
      • Riboldi M.
      • Orecchia R.
      • Baroni G.
      Real-time tumour tracking in particle therapy: technological developments and future perspectives.
      ,
      • Kubiak T.
      Particle therapy of moving targets-the strategies for tumour motion monitoring and moving targets irradiation.
      ,
      • Mori S.
      • Knopf A.C.
      • Umegaki K.
      Motion management in particle therapy.
      ]. 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. [
      • Phillips J.
      • Gueorguiev G.
      • Shackleford J.A.
      • Grassberger C.
      • Dowdell S.
      • Paganetti H.
      • et al.
      Computing proton dose to irregularly moving targets.
      ] 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 [
      • Meschini G.
      • Seregni M.
      • Molinelli S.
      • Vai A.
      • Phillips J.
      • Sharp G.C.
      • et al.
      Validation of a model for physical dose variations in irregularly moving targets treated with carbon ion beams.
      ,
      • Meschini G.
      • Kamp F.
      • Hofmaier J.
      • Reiner M.
      • Sharp G.
      • Paganetti H.
      • et al.
      Modeling RBE-weighted dose variations in irregularly moving abdominal targets treated with carbon ion 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 [
      • Meschini G.
      • Vai A.
      • Paganelli C.
      • Molinelli S.
      • Fontana G.
      • Pella A.
      • et al.
      Virtual 4DCT from 4DMRI for the management of respiratory motion in carbon ion therapy of abdominal tumors.
      ,
      • Bernatowicz K.
      • Peroni M.
      • Perrin R.
      • Weber D.C.
      • Lomax A.
      Four-Dimensional Dose Reconstruction for Scanned Proton Therapy Using Liver 4DCT-MRI.
      ]. Alternatively, image sequences can be generated with virtual phantoms, such as the 4D extended cardiac-torso (XCAT) phantom [
      • Segars W.P.
      • Sturgeon G.
      • Mendonca S.
      • Grimes J.
      • Tsui B.M.
      4D XCAT phantom for multimodality imaging research.
      ], 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.

      ,
      • Kostiukhina N.
      • Palmans H.
      • Stock M.
      • Georg D.
      • Knäusl B.
      Dynamic lung phantom commissioning for 4D dose assessment in proton therapy.
      ,
      • Kostiukhina N.
      • Palmans H.
      • Stock M.
      • Knopf A.
      • Georg D.
      • Knäusl B.
      Time-resolved dosimetry for validation of 4D dose calculation in PBS proton therapy.
      ]. 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 [
      • Paganetti H.
      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 [
      • Scholz M.
      • Kellerer A.M.
      • Kraft-Weyrather W.
      • Kraft G.
      Computation of cell survival in heavy ion beams for therapy: the model and its approximation.
      ,
      • Krämer M.
      • Scholz M.
      Rapid calculation of biological effects in ion radiotherapy.
      ,
      • Elsässer T.
      • Krämer M.
      • Scholz M.
      Accuracy of the local effect model for the prediction of biologic effects of carbon ion beams in vitro and in vivo.
      ,
      • Elsässer T.
      • Weyrather W.K.
      • Friedrich T.
      • Durante M.
      • Iancu G.
      • Krämer M.
      • et al.
      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 [
      • Scholz M.
      • Kellerer A.M.
      • Kraft-Weyrather W.
      • Kraft G.
      Computation of cell survival in heavy ion beams for therapy: the model and its approximation.
      ,
      • Krämer M.
      • Scholz M.
      Rapid calculation of biological effects in ion radiotherapy.
      ,
      • Krämer M.
      • Jäkel O.
      • Haberer T.
      • Kraft G.
      • Schardt D.
      • Weber U.
      Treatment planning for heavy-ion radiotherapy: Physical beam model and dose optimization.
      ,
      • Richter D.
      • Schwarzkopf A.
      • Trautmann J.
      • Krämer M.
      • Durante M.
      • Jäkel O.
      • et al.
      Upgrade and benchmarking of a 4D treatment planning system for scanned ion beam therapy.
      ,
      • Alliger C.
      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.

      2.1 The dose calculation algorithm

      The TRiP98 4D-dose calculation algorithm [
      • Richter D.
      • Schwarzkopf A.
      • Trautmann J.
      • Krämer M.
      • Durante M.
      • Jäkel O.
      • et al.
      Upgrade and benchmarking of a 4D treatment planning system for scanned ion beam therapy.
      ] 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) [
      • Chadwick K.H.
      • Leenhouts H.P.
      A molecular theory of cell survival.
      ,
      • McMahon S.J.
      The linear quadratic model: usage, interpretation and challenges.
      ]:
      S(D)=exp(-α¯·D-β¯·D2).
      (1)


      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 xi and yi in the iso-energetic slice (IES), which corresponds to the initial beam energy Bi. In order to take into account the entire particle spectrum, the dose averaged coefficients
      α¯=L-1lwldEdx(l)αl
      (2)


      β¯=L-1lwldEdx(l)βl
      (3)


      with the normalization factor
      L=lwldEdx(l)
      (4)


      are calculated based on the local effect model (LEM) [
      • Scholz M.
      • Kellerer A.M.
      • Kraft-Weyrather W.
      • Kraft G.
      Computation of cell survival in heavy ion beams for therapy: the model and its approximation.
      ] in the low dose approach. As described by Krämer and Scholz in [
      • Krämer M.
      • Scholz M.
      Rapid calculation of biological effects in ion radiotherapy.
      ], 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 wl represents the number of particles with these properties that contributes to the dose. The particle spectra in T and E and the resulting contributions Lik,αik and βik 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
      α¯=L-1iNspot,ikζikαik
      (5)


      β¯=L-1iNspot,ikζikβik
      (6)


      L=iNspot,ikζikLik.
      (7)


      Here, ζik,k=1,2 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
      DphysGy=C·LMeV/(g/cm)2
      (8)


      with C=1.602189·10-8. From Eq. 1 follows the biological effect
      -lnS=α¯Dphys+β¯Dphys2.
      (9)


      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:
      Dbio=-lnS/βX+(αX/(2βX))2-(αX/(2βX)).
      (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
      d(t)=Ak·sin2πTk·(t-t0,k)
      (11)


      with amplitude and period variations between subsequent periods, each period ranging from one minimum to another. Here, Ak and Tk denote the motion amplitude and period of the current motion and t0,k the time of the precedent minimum. The average period T was set to 5 s. Motion periods and amplitudes Ak were varied to uniformly distributed random values in the range [A-ΔA,A+ΔA] and [T-ΔT,T+ΔT], 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 5×106 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).
      Table 1Settings of the motion trajectory for the seven measurements.
      MotionA (mm)ΔA/A (%)T(s)ΔT/T (%)
      I10050
      II20050
      III30050
      IV202550
      V205050
      VI200550
      VII2050550
      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 [
      • Lis M.
      • Donetti M.
      • Newhauser W.
      • Durante M.
      • Dey J.
      • Weber U.
      • et al.
      A modular dose delivery system for treating moving targets with scanned ion beams: Performance and safety characteristics, and preliminary tests.
      ] 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 [
      • Aricò G.
      • Gehrke T.
      • Jakubek J.
      • Gallas R.
      • Berke S.
      • Jäkel O.
      • et al.
      Investigation of mixed ion fields in the forward direction for 220.5 MeV/u helium ion beams: Comparison between water and PMMA targets.
      ]) and the reference point of the PTW Octavius 1500XDR.
      A generalized 3D gamma analysis [
      • Wilson L.J.
      • Newhauser W.D.
      • Schneider C.W.
      An objective method to evaluate radiation dose distributions varying by three orders of magnitude.
      ] 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 Dthresh=0.15Gy, 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 ΔDA=0.08Gy 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 γ(3%/3mm) 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 N=25 virtual CTs was created using the software XCAT [
      • Segars W.P.
      • Sturgeon G.
      • Mendonca S.
      • Grimes J.
      • Tsui B.M.
      4D XCAT phantom for multimodality imaging research.
      ,
      • Segars W.P.
      • Mahesh M.
      • Beck T.J.
      • Frey E.C.
      • Tsui B.M.
      Realistic CT simulation using the 4D XCAT phantom.
      ]. 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. [
      • Lujan A.E.
      • Balter J.M.
      • Ten Haken R.K.
      A method for incorporating organ motion due to breathing into 3D dose calculations in the liver: Sensitivity to variations in motion.
      ] was selected. During breathing cycle k, so for t[t0,k;t0,k+Tk), the displacement d of the diaphragm’s dome from its standard end-exhale position in cranial direction was given by the equation
      d(t)=Ak·sin4πTk·(t-t0,k)+d0-vd·t.
      (12)


      Figure thumbnail gr2
      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, Ak and Tk represent motion amplitude of the diaphragm’s dome and period for respiratory cycle k, respectively. d0 is the start position of the first breathing cycle and vd is the velocity of a linear baseline drift. Three motion scenarios were considered for the dose calculation:
      • Periodic: Periodic motion from Eq. 12 with Ak=20mm and Tk=5s for all periods and no baseline drifts.
      • Fixed baseline: Nonperiodic motion, with Ak and Tk randomly generated with a uniform distribution between +-25% around A20mm and T=5s, and no baseline drift.
      • Baseline drift: The same motion as in the fixed baseline scenario, but with an additional linear baseline drift of vd=0.015mm/s.
      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 d0=12.5mm 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 d(t=0s)=d0, corresponding to CT number 7. For the dose calculation, this motion state was defined the reference state iref. The mean diaphragm displacement for each motion state j was calculated and the precomputed state i=i(j) 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 irefi for i>iref were directly exported from XCAT. For each voxel in the reference state iref, they define a vector to the corresponding position in the other states i. For i<iref, the transformation was obtained by mapping from state iref to the repetition of state i in a second, identical breathing cycle, i=i+N. 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 iref=iref+N was the end state. Now the vector field iiref 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 ji=i(j) 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 7,8,11,15,19,19,15,11,8 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 [
      • Ji Saitoh
      • Shirai K.
      • Mizukami T.
      • Abe T.
      • Ebara T.
      • Ohno T.
      • et al.
      Hypofractionated carbon-ion radiotherapy for stage I peripheral nonsmall cell lung cancer (GUNMA0701): Prospective phase II study.
      ] and hepatocellular carcinoma [
      • Kasuya G.
      • Kato H.
      • Yasuda S.
      • Tsuji H.
      • Yamada S.
      • Haruyama Y.
      • et al.
      Progressive hypofractionated carbon-ion radiotherapy for hepatocellular carcinoma: Combined analyses of 2 prospective trials.
      ]. The resulting dose distributions were compared by calculating the DVH metrics such as V95,D5 and D95 on the gross tumor volume. V95 is the fraction of the volume which receives a dose 95% of the prescribed dose and quantifies target coverage. D5 and D95 denote the dose that is exceeded in 5% and 95% of the target volume, respectively. The difference D5-D95 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.
      MotionAΔA/AΔT/Tγ pass rate (%)RtotRmedian
      (mm)(%)(%)3%/3 mm2%/2 mm
      I1000100.0088.140.99240.9553
      II2000100.0090.000.99130.946
      III3000100.0090.650.98890.93855
      IV20250100.0089.680.99150.94765
      V2050099.2088.000.98720.94805
      VI20050100.0084.130.98890.9431
      VII205050100.0086.780.99120.9528
      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 R2=0.99269 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 Rtot=0.9912. 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.
      Figure thumbnail gr3
      Fig. 3Recorded motion of the ionization chamber array detector for measurement VII. Black dots represent measured positions, the red line is an interpolation.
      Figure thumbnail gr4
      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.
      Figure thumbnail gr5
      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.
      Figure thumbnail gr6
      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 Rmedian 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.
      Figure thumbnail gr7
      Fig. 7Evolution of the correlation between the measured and reconstructed integral dose distributions over the course of the measurements listed in .
      Figure thumbnail gr8
      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 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.
      Figure thumbnail gr9
      Fig. 9Dependence of the γ(3%/3mm) pass rate on the duration of each reconstructed motion state.

      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 D5,D95 and V95 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 V95=83.7% (see Fig. 10f).
      Figure thumbnail gr10
      Fig. 10Calculated RBE-weighted dose distributions for the three motion scenarios described in (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.
      Table 3Results of the dose calculation for the three motion scenarios and two treatment modalities, where HI=D5-D95 depicts the homogeneity index.
      Motion scenario# RescansD5(%)D95(%)HI(%)V95(%)
      Periodic1106.990.016.981.2
      20101.696.94.799.6
      Fixed baseline1102.589.013.676.6
      20101.496.05.498.5
      Baseline drift1106.489.417.079.0
      20108.189.019.183.7

      4. Discussion

      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. [
      • Lis M.
      • Donetti M.
      • Newhauser W.
      • Durante M.
      • Dey J.
      • Weber U.
      • et al.
      A modular dose delivery system for treating moving targets with scanned ion beams: Performance and safety characteristics, and preliminary tests.
      ] 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 γ(3%/3mm) 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. [
      • Bert C.
      • Richter D.
      • Durante M.
      • Rietzel E.
      Scanned carbon beam irradiation of moving films: Comparison of measured and calculated response.
      ] 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 [
      • Kalantzopoulos C.
      • Meschini G.
      • Paganelli C.
      • Fontana G.
      • Vai A.
      • Preda L.
      • et al.
      Organ motion quantification and margins evaluation in carbon ion therapy of abdominal lesions.
      ].
      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 [
      • Bert C.
      • Richter D.
      • Durante M.
      • Rietzel E.
      Scanned carbon beam irradiation of moving films: Comparison of measured and calculated response.
      ], ICs in a 3D phantom [
      • Richter D.
      • Schwarzkopf A.
      • Trautmann J.
      • Krämer M.
      • Durante M.
      • Jäkel O.
      • et al.
      Upgrade and benchmarking of a 4D treatment planning system for scanned ion beam therapy.
      ] and cell survival studies [
      • Gemmel A.
      • Rietzel E.
      • Kraft G.
      • Durante M.
      • Bert C.
      Calculation and experimental verification of the RBE-weighted dose for scanned ion beams in the presence of target motion.
      ]. 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.

      ,
      • Kostiukhina N.
      • Palmans H.
      • Stock M.
      • Georg D.
      • Knäusl B.
      Dynamic lung phantom commissioning for 4D dose assessment in proton therapy.
      ,
      • Kostiukhina N.
      • Palmans H.
      • Stock M.
      • Knopf A.
      • Georg D.
      • Knäusl B.
      Time-resolved dosimetry for validation of 4D dose calculation in PBS proton therapy.
      ] 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 [
      • Bernatowicz K.
      • Peroni M.
      • Perrin R.
      • Weber D.C.
      • Lomax A.
      Four-Dimensional Dose Reconstruction for Scanned Proton Therapy Using Liver 4DCT-MRI.
      ]. The feasibility to create suitable data sets was demonstrated for another dose calculation algorithm for carbon ion beams [
      • Meschini G.
      • Kamp F.
      • Hofmaier J.
      • Reiner M.
      • Sharp G.
      • Paganetti H.
      • et al.
      Modeling RBE-weighted dose variations in irregularly moving abdominal targets treated with carbon ion beams.
      ,
      • Meschini G.
      • Vai A.
      • Paganelli C.
      • Molinelli S.
      • Fontana G.
      • Pella A.
      • et al.
      Virtual 4DCT from 4DMRI for the management of respiratory motion in carbon ion therapy of abdominal tumors.
      ], 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 [
      • Eiben B.
      • Bertholet J.
      • Menten M.J.
      • Nill S.
      • Oelfke U.
      • McClelland J.R.
      Consistent and invertible deformation vector fields for a breathing anthropomorphic phantom: A post-processing framework for the XCAT phantom.
      ] 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 [
      • Graeff C.
      • Constantinescu A.
      • Lüchtenborg R.
      • Durante M.
      • Bert C.
      Multigating, a 4D optimized beam tracking in scanned ion beam therapy.
      ,
      • Wolf M.
      • Anderle K.
      • Durante M.
      • Graeff C.
      Robust treatment planning with 4D intensity modulated carbon ion therapy for multiple targets in stage IV non-small cell lung cancer.
      ,
      • Graeff C.
      Motion mitigation in scanned ion beam therapy through 4D-optimization.
      ] was proven to deliver homogeneous and more conformal doses faster for moving targets [
      • Lis M.
      • Donetti M.
      • Newhauser W.
      • Durante M.
      • Dey J.
      • Weber U.
      • et al.
      A modular dose delivery system for treating moving targets with scanned ion beams: Performance and safety characteristics, and preliminary tests.
      ]. Developments to combine this technique with residual real time tracking are ongoing [

      LeDoeuff G. Beam tracking for compensation of irregular motion (MSc Thesis). Darmstadt 2018.

      ]. 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. [
      • Takao S.
      • Miyamoto N.
      • Matsuura T.
      • Onimaru R.
      • Katoh N.
      • Inoue T.
      • 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.
      ]. Also sudden changes of the baseline position were observed in real patient motion [
      • Kalantzopoulos C.
      • Meschini G.
      • Paganelli C.
      • Fontana G.
      • Vai A.
      • Preda L.
      • et al.
      Organ motion quantification and margins evaluation in carbon ion therapy of abdominal lesions.
      ]. 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 [
      • Riboldi M.
      • Orecchia R.
      • Baroni G.
      Real-time tumour tracking in particle therapy: technological developments and future perspectives.
      ,
      • Bert C.
      • Herfarth K.
      Management of organ motion in scanned ion beam therapy.
      ], 4D optimization [
      • Graeff C.
      • Constantinescu A.
      • Lüchtenborg R.
      • Durante M.
      • Bert C.
      Multigating, a 4D optimized beam tracking in scanned ion beam therapy.
      ,
      • Wolf M.
      • Anderle K.
      • Durante M.
      • Graeff C.
      Robust treatment planning with 4D intensity modulated carbon ion therapy for multiple targets in stage IV non-small cell lung cancer.
      ,
      • Graeff C.
      Motion mitigation in scanned ion beam therapy through 4D-optimization.
      ] and a combination of both [

      LeDoeuff G. Beam tracking for compensation of irregular motion (MSc Thesis). Darmstadt 2018.

      ] 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.

      References

        • Durante M.
        • Loeffler J.S.
        Charged particles in radiation oncology.
        Nat Rev Clin Oncol. 2010; 7: 37-43https://doi.org/10.1038/nrclinonc.2009.183
        • Grau C.
        • Durante M.
        • Georg D.
        • Langendijk J.A.
        • Weber D.C.
        Particle therapy in Europe.
        Mol Oncol. 2020; 14: 1492-1499https://doi.org/10.1002/1878-0261.12677
        • Riboldi M.
        • Orecchia R.
        • Baroni G.
        Real-time tumour tracking in particle therapy: technological developments and future perspectives.
        Lancet Oncol. 2012; 13https://doi.org/10.1016/S1470-2045(12)70243-7
        • Bert C.
        • Durante M.
        Motion in radiotherapy: particle therapy.
        Phys Med Biol. 2011; 56https://doi.org/10.1088/0031-9155/56/16/R01
        • He P.
        • Mori S.
        Perturbation analysis of 4d dose distribution for scanned carbon-ion beam radiotherapy.
        Phys Med. 2020; 74: 74-82https://doi.org/10.1016/j.ejmp.2020.05.003
        • Takao S.
        • Miyamoto N.
        • Matsuura T.
        • Onimaru R.
        • Katoh N.
        • Inoue T.
        • 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.
        Int J Radiat Oncol Biol Phys. 2016; 94: 172-180https://doi.org/10.1016/j.ijrobp.2015.09.024
        • Dhont J.
        • Vandemeulebroucke J.
        • Burghelea M.
        • Poels K.
        • Depuydt T.
        • Van Den Begin R.
        • et al.
        The long- and short-term variability of breathing induced tumor motion in lung and liver over the course of a radiotherapy treatment.
        Radiother Oncol. 2018; 126: 339-346https://doi.org/10.1016/j.radonc.2017.09.001
        • Bert C.
        • Herfarth K.
        Management of organ motion in scanned ion beam therapy.
        Radiat Oncol. 2017; 12: 1-6https://doi.org/10.1186/s13014-017-0911-z
        • Graeff C.
        • Constantinescu A.
        • Lüchtenborg R.
        • Durante M.
        • Bert C.
        Multigating, a 4D optimized beam tracking in scanned ion beam therapy.
        Technol Cancer Res Treat. 2014; 13: 497-504https://doi.org/10.7785/tcrtexpress.2013.600277
      1. 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.

        • Wolf M.
        • Anderle K.
        • Durante M.
        • Graeff C.
        Robust treatment planning with 4D intensity modulated carbon ion therapy for multiple targets in stage IV non-small cell lung cancer.
        Phys Med Biol. 2020; 65215012https://doi.org/10.1088/1361-6560/aba1a3
        • Kanai T.
        • Paz A.
        • Furuichi W.
        • Liu C.S.
        • He P.
        • Mori S.
        Four-dimensional carbon-ion pencil beam treatment planning comparison between robust optimization and range-adapted internal target volume for respiratory-gated liver and lung treatment.
        Phys Med. 2020; 80https://doi.org/10.1016/j.ejmp.2020.11.009
        • Meijers A.
        • Jakobi A.
        • Stützer K.
        • Guterres Marmitt G.
        • Both S.
        • Langendijk J.A.
        • et al.
        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.
        Med Phys. 2019; 46: 1140-1149https://doi.org/10.1002/mp.13371
        • Vedam S.S.
        • Keall P.J.
        • Kini V.R.
        • Mostafavi H.
        • Shukla H.P.
        • Mohan R.
        Acquiring a four-dimensional computed tomography dataset using an external respiratory signal.
        Phys Med Biol. 2003; 48: 45-62https://doi.org/10.1088/0031-9155/48/1/304
        • Meijers A.
        • Knopf A.C.
        • Crijns A.P.
        • Ubbels J.F.
        • Niezink A.G.
        • Langendijk J.A.
        • et al.
        Evaluation of interplay and organ motion effects by means of 4D dose reconstruction and accumulation.
        Radiother Oncol. 2020; 150: 268-274https://doi.org/10.1016/j.radonc.2020.07.055
        • Li R.
        • Jia X.
        • Lewis J.H.
        • Gu X.
        • Folkerts M.
        • Men C.
        • et al.
        Real-time volumetric image reconstruction and 3D tumor localization based on a single x-ray projection image for lung cancer radiotherapy.
        Med Phys. 2010; 37: 2822-2826https://doi.org/10.1118/1.3426002
        • Raaymakers B.W.
        • Raaijmakers A.J.
        • Lagendijk J.J.
        Feasibility of MRI guided proton therapy: magnetic field dose effects.
        Phys Med Biol. 2008; 53: 5615-5622https://doi.org/10.1088/0031-9155/53/20/003
        • O’Shea T.
        • Bamber J.
        • Fontanarosa D.
        • Van Der Meer S.
        • Verhaegen F.
        • Harris E.
        Review of ultrasound image guidance in external beam radiotherapy part II: Intra-fraction motion management and novel applications.
        Phys Med Biol. 2016; 61: R90-R137https://doi.org/10.1088/0031-9155/61/8/R90
        • Berbeco R.I.
        • Jiang S.B.
        • Sharp G.C.
        • Chen G.T.
        • Mostafavi H.
        • Shirato H.
        Integrated radiotherapy imaging system (IRIS): design considerations of tumour tracking with linac gantry-mounted diagnostic x-ray systems with flat-panel detectors.
        Phys Med Biol. 2004; 49: 243-255https://doi.org/10.1088/0031-9155/49/2/005
        • Dong B.
        • Graves Y.J.
        • Jia X.
        • Jiang S.B.
        Optimal surface marker locations for tumor motion estimation in lung cancer radiotherapy.
        Phys Med Biol. 2012; 57: 8201-8215https://doi.org/10.1088/0031-9155/57/24/8201
        • Kubiak T.
        Particle therapy of moving targets-the strategies for tumour motion monitoring and moving targets irradiation.
        Br J Radiol. 2016; 89https://doi.org/10.1259/bjr.20150275
        • Mori S.
        • Knopf A.C.
        • Umegaki K.
        Motion management in particle therapy.
        Med Phys. 2018; 45: e994-e1010https://doi.org/10.1002/mp.12679
        • Phillips J.
        • Gueorguiev G.
        • Shackleford J.A.
        • Grassberger C.
        • Dowdell S.
        • Paganetti H.
        • et al.
        Computing proton dose to irregularly moving targets.
        Phys Med Biol. 2014; 59https://doi.org/10.1088/0031-9155/59/15/4261
        • Meschini G.
        • Seregni M.
        • Molinelli S.
        • Vai A.
        • Phillips J.
        • Sharp G.C.
        • et al.
        Validation of a model for physical dose variations in irregularly moving targets treated with carbon ion beams.
        Med Phys. 2019; 46: 3663-3673https://doi.org/10.1002/mp.13662
        • Meschini G.
        • Kamp F.
        • Hofmaier J.
        • Reiner M.
        • Sharp G.
        • Paganetti H.
        • et al.
        Modeling RBE-weighted dose variations in irregularly moving abdominal targets treated with carbon ion beams.
        Med Phys. 2020; 47: 2768-2778https://doi.org/10.1002/mp.14135
        • Meschini G.
        • Vai A.
        • Paganelli C.
        • Molinelli S.
        • Fontana G.
        • Pella A.
        • et al.
        Virtual 4DCT from 4DMRI for the management of respiratory motion in carbon ion therapy of abdominal tumors.
        Med Phys. 2020; 47: 909-916https://doi.org/10.1002/mp.13992
        • Bernatowicz K.
        • Peroni M.
        • Perrin R.
        • Weber D.C.
        • Lomax A.
        Four-Dimensional Dose Reconstruction for Scanned Proton Therapy Using Liver 4DCT-MRI.
        Int J Radiat Oncol Biol Phys. 2016; 95: 216-223https://doi.org/10.1016/j.ijrobp.2016.02.050
        • Segars W.P.
        • Sturgeon G.
        • Mendonca S.
        • Grimes J.
        • Tsui B.M.
        4D XCAT phantom for multimodality imaging research.
        Med Phys. 2010; 37: 4902-4915https://doi.org/10.1118/1.3480985
      2. 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.

        • Kostiukhina N.
        • Palmans H.
        • Stock M.
        • Georg D.
        • Knäusl B.
        Dynamic lung phantom commissioning for 4D dose assessment in proton therapy.
        Phys Med Biol. 2019; 64235001https://doi.org/10.1088/1361-6560/ab5132
        • Kostiukhina N.
        • Palmans H.
        • Stock M.
        • Knopf A.
        • Georg D.
        • Knäusl B.
        Time-resolved dosimetry for validation of 4D dose calculation in PBS proton therapy.
        Phys Med Biol. 2020; 65125015https://doi.org/10.1088/1361-6560/ab8d79
        • Paganetti H.
        Relative biological effectiveness (RBE) values for proton beam therapy. Variations as a function of biological endpoint, dose, and linear energy transfer.
        Phys Med Biol. 2014; 59https://doi.org/10.1088/0031-9155/59/22/R419
        • Scholz M.
        • Kellerer A.M.
        • Kraft-Weyrather W.
        • Kraft G.
        Computation of cell survival in heavy ion beams for therapy: the model and its approximation.
        Radiat Environ Biophys. 1997; 36: 59-66https://doi.org/10.1007/s004110050055
        • Krämer M.
        • Scholz M.
        Rapid calculation of biological effects in ion radiotherapy.
        Phys Med Biol. 2006; 51: 1959-1970https://doi.org/10.1088/0031-9155/51/8/001
        • Elsässer T.
        • Krämer M.
        • Scholz M.
        Accuracy of the local effect model for the prediction of biologic effects of carbon ion beams in vitro and in vivo.
        Int J Radiat Oncol Biol Phys. 2008; 71https://doi.org/10.1016/j.ijrobp.2008.02.037
        • Elsässer T.
        • Weyrather W.K.
        • Friedrich T.
        • Durante M.
        • Iancu G.
        • Krämer M.
        • et al.
        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.
        Int J Radiat Oncol Biol Phys. 2010; 78https://doi.org/10.1016/j.ijrobp.2010.05.014
        • Krämer M.
        • Jäkel O.
        • Haberer T.
        • Kraft G.
        • Schardt D.
        • Weber U.
        Treatment planning for heavy-ion radiotherapy: Physical beam model and dose optimization.
        Phys Med Biol. 2000; 45: 3299-3317https://doi.org/10.1088/0031-9155/45/11/313
        • Richter D.
        • Schwarzkopf A.
        • Trautmann J.
        • Krämer M.
        • Durante M.
        • Jäkel O.
        • et al.
        Upgrade and benchmarking of a 4D treatment planning system for scanned ion beam therapy.
        Med Phys. 2013; 40https://doi.org/10.1118/1.4800802
        • Alliger C.
        Zeitlich aufgelöste Dosisberechnung für gescannte Kohlenstoffionen auf sequentiellen Bilddaten zur Erfassung unregelmäßiger Bewegungen (Master Thesis).
        Darmstadt. 2018;
        • Chadwick K.H.
        • Leenhouts H.P.
        A molecular theory of cell survival.
        Phys Med Biol. 1973; 18: 78-87https://doi.org/10.1088/0031-9155/18/1/007
        • McMahon S.J.
        The linear quadratic model: usage, interpretation and challenges.
        Phys Med Biol. 2019; 64https://doi.org/10.1088/1361-6560/aaf26a
        • Lis M.
        • Donetti M.
        • Newhauser W.
        • Durante M.
        • Dey J.
        • Weber U.
        • et al.
        A modular dose delivery system for treating moving targets with scanned ion beams: Performance and safety characteristics, and preliminary tests.
        Phys Med. 2020; 76: 307-316https://doi.org/10.1016/j.ejmp.2020.07.029
        • Aricò G.
        • Gehrke T.
        • Jakubek J.
        • Gallas R.
        • Berke S.
        • Jäkel O.
        • et al.
        Investigation of mixed ion fields in the forward direction for 220.5 MeV/u helium ion beams: Comparison between water and PMMA targets.
        Phys Med Biol. 2017; 62https://doi.org/10.1088/1361-6560/aa875e
        • Wilson L.J.
        • Newhauser W.D.
        • Schneider C.W.
        An objective method to evaluate radiation dose distributions varying by three orders of magnitude.
        Med Phys. 2019; 46: 1888-1895https://doi.org/10.1002/mp.13420
        • Segars W.P.
        • Mahesh M.
        • Beck T.J.
        • Frey E.C.
        • Tsui B.M.
        Realistic CT simulation using the 4D XCAT phantom.
        Med Phys. 2008; 35: 3800-3808https://doi.org/10.1118/1.2955743
        • Lujan A.E.
        • Balter J.M.
        • Ten Haken R.K.
        A method for incorporating organ motion due to breathing into 3D dose calculations in the liver: Sensitivity to variations in motion.
        Med Phys. 2003; 30: 2643-2649https://doi.org/10.1118/1.1609057
        • Ji Saitoh
        • Shirai K.
        • Mizukami T.
        • Abe T.
        • Ebara T.
        • Ohno T.
        • et al.
        Hypofractionated carbon-ion radiotherapy for stage I peripheral nonsmall cell lung cancer (GUNMA0701): Prospective phase II study.
        Cancer Med. 2019; 8https://doi.org/10.1002/cam4.2561
        • Kasuya G.
        • Kato H.
        • Yasuda S.
        • Tsuji H.
        • Yamada S.
        • Haruyama Y.
        • et al.
        Progressive hypofractionated carbon-ion radiotherapy for hepatocellular carcinoma: Combined analyses of 2 prospective trials.
        Cancer. 2017; 123https://doi.org/10.1002/cncr.30816
        • Bert C.
        • Richter D.
        • Durante M.
        • Rietzel E.
        Scanned carbon beam irradiation of moving films: Comparison of measured and calculated response.
        Radiat Oncol. 2012; 7https://doi.org/10.1186/1748-717X-7-55
        • Kalantzopoulos C.
        • Meschini G.
        • Paganelli C.
        • Fontana G.
        • Vai A.
        • Preda L.
        • et al.
        Organ motion quantification and margins evaluation in carbon ion therapy of abdominal lesions.
        Phys Med. 2020; 75: 33-39https://doi.org/10.1016/j.ejmp.2020.05.014
        • Gemmel A.
        • Rietzel E.
        • Kraft G.
        • Durante M.
        • Bert C.
        Calculation and experimental verification of the RBE-weighted dose for scanned ion beams in the presence of target motion.
        Phys Med Biol. 2011; 56: 7337-7351https://doi.org/10.1088/0031-9155/56/23/001
        • Eiben B.
        • Bertholet J.
        • Menten M.J.
        • Nill S.
        • Oelfke U.
        • McClelland J.R.
        Consistent and invertible deformation vector fields for a breathing anthropomorphic phantom: A post-processing framework for the XCAT phantom.
        Phys Med Biol. 2020; 65https://doi.org/10.1088/1361-6560/ab8533
        • Graeff C.
        Motion mitigation in scanned ion beam therapy through 4D-optimization.
        Phys Med. 2014; 30: 570-577https://doi.org/10.1016/j.ejmp.2014.03.011
      3. LeDoeuff G. Beam tracking for compensation of irregular motion (MSc Thesis). Darmstadt 2018.