Deformable image registration of the treatment planning CT with proton radiographies in perspective of adaptive proton therapy

The purpose of this work is to investigate the potentiality of using a limited number of in-room proton radiographies to compensate anatomical changes in adaptive proton therapy. The treatment planning CT is adapted to the treatment delivery scenario relying on 2D-3D deformable image registration (DIR). The proton radiographies, expressed in water equivalent thickness (WET) are simulated for both list-mode and integration-mode detector configurations in pencil beam scanning. Geometrical and analytical simulations of an anthropomorphic phantom in the presence of anatomical changes due to breathing are adopted. A Monte Carlo simulation of proton radiographies based on a clinical CT image in the presence of artificial anatomical changes is also considered. The accuracy of the 2D-3D DIR, calculated as root mean square error, strongly depends on the considered anatomical changes and is considered adequate for promising adaptive proton therapy when comparable to the accuracy of conventional 3D-3D DIR. In geometrical simulation, this is achieved with a minimum of eight/nine radiographies (more than 90% accuracy). Negligible improvement (sim1%) is obtained with the use of 180 radiographies. Comparing different detector configurations, superior accuracy is obtained with list-mode than integration-mode max (WET with maximum occurrence) and mean (average WET weighted by occurrences). Moreover, integration-mode max performs better than integration-mode mean. Results are minimally affected by proton statistics. In analytical simulation, the anatomical changes are approximately compensated (about 60%–70% accuracy) with two proton radiographies and minor improvement is observed with nine proton radiographies. In clinical data, two proton radiographies from list-mode have demonstrated better performance than nine from integration-mode (more than 100% and about 50%–70% accuracy, respectively), even avoiding the finer grid spacing of the last numerical optimization stage. In conclusion, the choice of detector configuration as well as the amount and complexity of the considered anatomical changes determine the minimum number of radiographies to be used.


Introduction
Prior to treatment delivery, the patient position in the treatment room is typically corrected based on rigid registration of the treatment planning CT to the in-room x-ray radiographic imaging. The registration is implemented as a numerical optimization of the matching between the forward-projection of the treatment planning CT to the radiographies (Fattori et al 2015). The forward-projection is calculated as integration of the treatment planning CT along the photon beam direction, the so-called digitally reconstructed radiograph (DRR). The approach is referred to as two-dimensional to three-dimensional (2D-3D) registration, where the moving image is the treatment planning CT, whereas the fixed images are the x-ray radiographies acquired in the treatment delivery scenario. Relying on a minimum of two radiographies, a six-parameter correction vector is obtained. The patient position is then corrected according to the inverse vector.
Relying on deformable image registration (DIR), anatomical changes can be accounted for by adapting the treatment planning CT to the treatment delivery scenario, according to the concept of adaptive radiation therapy (ART) (Yan et al 1997). When an in-room tomographic imaging is available, conventional 3D-3D DIR can be applied (Peroni et al 2012, Landry et al 2015. In photon beam therapy, commercial cone beam CT (CBCT) devices are nowadays commonly integrated in the treatment room by mounting the x-ray imaging systems on the rotating photon beam gantry, perpendicular to the therapeutic beam direction (Létourneau et al 2005). In ion beam therapy, rotating gantries are mostly designed for proton beams, whereas integrated CBCT is just entering the clinical workflow (Mcdonough and Tinnel 2007, Dupont 2011, Landry et al 2018. However, a rotating proton beam gantry able to deliver highly energetic proton beams can be directly exploited for emerging ion imaging systems (Johnson et al 2017). By mounting a dedicated detector on the proton beam gantry, a proton tomography (pCT) can be in principle obtained by means of tomographic image reconstruction of multiple proton radiographies. pCT can represent therefore the native in-room tomographic imaging solution for ART in proton therapy. Proton imaging systems are differently conceived for detectors capable of measuring single protons (list-mode) or just multiple protons of single pencil beams (integration-mode). List-mode detectors are typically composed by two trackers (Bashkirov et al 2016, Taylor et al 2016, Mattiazzo et al 2018, upstream and downstream with respect to the patient, and an absorption detector to measure the residual energy of each proton. Integration-mode detectors are composed of only the absorption detector, which retrieves the mixed range of the pencil beam through the time-resolved energy loss at multiple beam energies (Testa et al 2013, Würl et al 2020 or through multiple detection layers (Rinaldi et al 2013, Meyer et al 2017. The measurement is then converted into water equivalent thickness (WET), which corresponds to the integral stopping power of the protons relative to water, or relative stopping power (RSP). Hence, the RSP of the patient anatomy can be reconstructed in pCT.
The possibility to perform pCT imaging is constrained by the clinical requirement of imaging dose minimization and the geometrical and technical limitations of integrating the detectors (especially for list-mode) in the treatment room. Moreover, tomographic image reconstruction of a limited number of proton radiographies compromises the accuracy of the 3D-3D DIR. Therefore, the 2D-3D DIR, as enabled by a limited number of proton radiographies, is worth investigating according to the perspective of ART in proton therapy. Similar to patient position correction, the 2D-3D DIR is implemented as a numerical optimization of the matching between the forward-projection of the treatment planning CT, calibrated to RSP, to the proton radiographies, expressing the WET. The forward-projection is calculated as integration of the calibrated treatment planning CT along the single proton or central pencil beam trajectory/axis. For list-mode radiographies, a proton trajectory is estimated, whereas for integration-mode radiographies the straight (most probable) central pencil beam axis is assumed. Relying on a minimum of two proton radiographies, a deformation field in image domain can be obtained (Long et al 2010, Brock et al 2010. The deformation field is then applied to the calibrated treatment planning CT, thus compensating for the anatomical changes. Alternatively, a deformation field in projection domain can be obtained based on 2D-2D DIR (Gianoli et al 2016a), and subsequently mapped in image domain (Gianoli et al 2014). However, this approach is affected by inaccuracies in domain conversions (Gianoli et al 2016b).
In this work, a recently proposed 2D-3D DIR (Palaniappan et al 2019) is considered for ART in proton therapy. The compensation of the anatomical changes between the treatment planning and treatment delivery scenario (figure 1) is extensively investigated in an anthropomorphic phantom. Relying on a purely geometrical simulation, the parameters of the 2D-3D DIR are optimized and the performance is assessed for different number of radiographies. The minimum number of radiographies able to accurately compensate the anatomical changes is identified. Subsequently, a recently proposed analytical simulation of proton imaging is adopted to compare the performance of the 2D-3D DIR when applied to list-mode and integration-mode proton radiographies relying on proton pencil beam scanning (Gianoli et al 2019) (figure 2(A)). The influence of proton statistics and the number of radiographies are investigated to define the potentialities (and limitations) of list-mode and integration-mode detector configurations. Finally, the 2D-3D DIR is preliminary applied to more realistic Monte Carlo simulation of proton radiographies based on clinical data (head and neck), for both list-mode and integration-mode detector configurations, to identify the challenges for future studies.

2D-3D DIR
A quasi-Newton Broyden-Fletcher-Goldfarb-Shanno optimization algorithm, commonly employed in 3D-3D DIR, is adopted. The optimization algorithm is based on Newton direction for line search calculating the objective function around the minimum as a quadratic Taylor polynomial and using an approximated Hessian matrix (hence, quasi-Newton) (Dennis Jr and Moré 1977). With this optimization algorithm, the deformation field is described by a free-form deformation based on cubic basis splines (B-splines) (Rueckert et al 1999, Kroon 2011. The basic idea of the free-form deformation B-spline parametrization is to deform the moving image by handling control points that produce a smooth and continuous deformation field in the image domain. The optimization algorithm, extended to 2D-3D DIR (Palaniappan et al 2019) (figures 1 and 2(B)), is applied to proton radiographies in geometrical and analytical simulations based on anthropomorphic phantom and Monte Carlo simulation based on clinical data. A stage of the optimization algorithm is defined by the spacing between control points, referred to as grid spacing, and by the number of iterations. More than one stage with the same number of iterations is considered and for subsequent stages, the spacing is halved, aiming at progressive refinement of the produced deformation field. The deformation field is obtained by matching the proton radiographies to the forward-projection of the treatment planning CT, calibrated to RSP. At each iteration, the deformation field is applied to the calibrated treatment planning CT, thus compensating for the anatomical changes. Hence, the calculation of the proton DRRs is embedded in the optimization algorithm ( figure 2(B)). The objective function in radiography domain is then calculated between the radiographies and the proton DRRs. The root mean square error (RMSE), and the normalized mutual information (NMI), which is expected to be less sensitive to possible inaccuracies of calibration curves, are adopted to quantify the matching between radiographies and proton DRRs.

Anthropomorphic phantom
The inhale and exhale breathing phases of an anthropomorphic phantom (Segars et al 2002) are considered to simulate large and complex inter-fractional anatomical changes between the ground truth pCT image and the treatment planning CT. The breathing of the anthropomorphic phantom is modeled by non-uniform rational basis splines relying on a superior-inferior diaphragm motion and anterior-posterior chest wall expansion curves. The image size of the phantom is 128 × 128 × 100 voxels and the voxel size is 3 × 3 × 3 mm 3 . First, the original photon linear attenuation coefficients of the inhale and exhale breathing phases are linearly converted to Hounsfield unit (HU). Subsequently, the inhale phase is converted to RSP, relying on the true calibration curve, thus generating the ground truth pCT (Gianoli et al 2019).

Geometrical simulation: ideal proton radiographies and proton DRRs of the anthropomorphic phantom
The ideal proton radiographies are modeled as the forward-projection of the ground truth pCT calculating the WET as the integral RSP along a straight line. The forward-projection is based on parallel straight lines spaced by the voxel size of the anthropomorphic phantom. The size of the radiography is therefore 128 × 100. The number of radiographies is set equal to 180, uniformly covering 180 • . The calculation of the ideal proton DRRs is based on the forward-projection, but applied to the calibrated treatment planning CT. The number of stages is set equal to three, according to 64, 32 and 16 voxels as grid spacing. The number of iterations for each stage is set equal to 15.

Analytical simulation: detector-dependent proton radiographies and proton DRRs of the anthropomorphic phantom
Individual proton trajectories in pencil beam scanning are analytically simulated and the proton radiographies are calculated as forward-projection of the ground truth pCT along these trajectories . The parallel straight central pencil beam axes are spaced by the voxel size of the anthropomorphic phantom. The curved proton trajectories are traced within the statistical distribution of a multiple Coulomb scattering model, originally given in uniform water and extended to non-uniform water equivalent materials, neglecting energy straggling and nuclear interactions. Therefore, the simulated proton trajectory deviates (due to multiple Coulomb scattering) from the central pencil beam axis in function of the RSP of the anthropomorphic phantom. The WET is therefore calculated as the integral RSP along the simulated proton trajectory. Each pencil beam with an energy of 280 MeV is composed of 25, 50, 75 and 100 protons per pencil beam (referred to as proton statistics). The 180 • are uniformly covered by nine radiographies plus the lateral radiography (at zero degrees) that complement the frontal one (at 90 • ), thus maximizing the orthogonality between radiographies. The proton radiographies are generated according to ideal list-mode and integration-mode detector configurations (figure 2(A)), implying no limitation in temporal resolution, WET resolution of the absorption detector as well as spatial resolution of the tracker. The list-mode proton radiography is composed by the WET for each proton of the pencil beam. For integration-mode proton radiographies, the exact solution of the mixed range signal, requiring linear decomposition in realistic integration-mode detectors (Meyer et al 2017), is considered. For this purpose, the WET of each proton is assumed to contribute to the WET histogram for each pencil beam. The WET histogram for each pencil beam offers a weighted mean WET (individual WET values weighted by the occurrences) and a maximum WET (the WET value with the maximum occurrence), thus referring to integration-mode mean and max radiographies, respectively. The proton DRRs are calculated for both list-mode and integration-mode detector configurations. For list-mode radiographies, the proton DRRs are calculated as integration of the calibrated treatment planning CT along the estimated proton trajectories (proton-wise). The proton trajectories are estimated based on the scattering model in uniform water, as the RSP of the anthropomorphic phantom is actually unknown. The estimation of the proton trajectories is relevant to the proton radiographies and thus, to the treatment delivery scenario. Therefore, the proton DRRs for list-mode detector configuration are affected by inaccuracies due to the anatomical changes. For integration-mode mean and max, the proton DRRs are calculated as integration of the calibrated treatment planning CT along the straight central pencil beam axis (pencil beam-wise). The central pencil beam axis is therefore assumed as the estimated proton trajectory for both integration-mode mean and max. The schematic representation of the proposed 2D-3D DIR, with emphasis on detector-dependent proton radiographies and proton DRRs, is shown in figure 2(B). The number of stages is set equal to two, according to 64 and 32 voxels as grid spacing. The number of iterations for each stage is set equal to 15.

Monte Carlo simulation of detector-dependent proton radiographies and proton DRRs of clinical data
A clinical CT image of a head and neck patient from the Department of Radiation Oncology at the Universitätsklinikum der Ludwig-Maximilians-Universität München is considered. The CT is converted to the ground truth pCT relying on a clinical monotonic calibration curve . The CT image size is 314 × 314 × 10 voxels, the pixel size is 0.1074 × 0.1074 cm 2 and the slice thickness is 0.3 cm. The ground truth pCT is adopted in Monte Carlo simulation of proton radiographies from ideal list-mode and integration-mode detector configurations in pencil beam scanning. The proton statistics is set equal to 400 protons per pencil beam (energy 199.94 MeV, pencil beam size 8.5 mm). The proton interaction and transport are simulated in FLUKA (Ferrari et al 2005, Böhlen et al 2014, relying on a customized simulation framework . The number of radiographies is set equal to nine, uniformly covering 180 • . List-mode and integration-mode proton radiographies and proton DRRs are obtained as for the analytical simulation ( figure 2(A)).
Artificial anatomical changes are obtained by superimposing two Gaussians of opposite signs for each component of the deformation field. The deformation fields are laterally placed on different positions. The lateral component of the deformation field has therefore opposite signs to simulate weight loss. The Gaussian amplitude and standard deviation are ±8 and 16 mm, respectively. The deformation field is then applied to the ground truth pCT to obtain the calibrated treatment planning CT. For integration-mode radiographies the stages are set equal to three (64, 32 and 16 grid spacing). Because of computational limitations, the number of stages is instead reduced to two (64 and 32 as grid spacing) for list-mode radiographies.

Performance quantification
The performance of the 2D-3D DIR is evaluated in comparison to the conventional 3D-3D DIR of the calibrated treatment planning CT to the ground truth pCT, relying on the same registration parameters (grid spacing and number of iterations). The objective function in radiography domain for the conventional 3D-3D DIR is the RMSE.
Relying on the ground truth pCT, the RMSE of the RSP is adopted as metric for accuracy quantification in image domain (intensity-based). For the anthropomorphic phantom, a distance-based metric (percentage overlap of anatomical regions of interest), prior to and after 2D-3D DIR, is also considered. The chosen anatomical regions of interest are lung lesion, lung, liver and heart. The deformation field is applied to the

Geometrical simulation of the anthropomorphic phantom
The comparison between the RMSE and the NMI as objective functions is performed on the geometrical simulation and reported in figure 3 for different number of proton radiographies uniformly covering the 180 • . Relying on the RMSE of RSP (table 1), at least eight or nine proton radiographies are needed to approach the conventional 3D-3D DIR and to match (~1%) the 2D-3D DIR based on 180 proton radiographies. The Dice similarity coefficients for each segmented anatomical region of interest prior to and after 2D-3D DIR in comparison to the conventional 3D-3D DIR at the end of the last stages are reported in table 2. The Dice similarity coefficients for each stage during the 2D-3D DIR and 3D-3D DIR are shown in the appendix ( figure A1). Compared to the RSME of RSP, the Dice similarity coefficients demonstrate less sensitivity to the number of proton radiographies. However, the Dice similarity coefficients strongly depend on the size, shape (simple/complex) and position with respect to the diaphragm of the region of interest. The latter correlates with the amplitude of the applied anatomical changes.

Analytical simulation of the anthropomorphic phantom
The comparison between different proton statistics is performed for analytical simulation of integration-mode mean (figure 4), integration-mode max (figure 5) and list-mode (figure 6) for different number of proton radiographies. In figure 7, the overlay of the ground truth pCT and the calibrated treatment planning CT, prior and after 2D-3D DIR, is shown for analytical simulation of anthropomorphic phantom. The Dice similarity coefficients for segmented anatomical regions of interest are reported in the appendix (figures A2, A3 and A4). The superior performance of the NMI over the RMSE as the objective function is confirmed. More stable results are obtained with NMI rather than RMSE as the objective function, especially appreciable in integration-mode max. The sharper interfaces of the proton radiographies generated by the selection of the WET with the maximum occurrence strongly affect the RMSE as the   objective function. This selection is prone to inaccuracies as proton statistics decrease. However, this feature in combination with the consideration of only two proton radiographies enables the 2D-3D DIR to broaden the searching space of the 2D-3D DIR, thus resulting in slightly better accuracy for the lower proton statistics than the higher ones for both RMSE and NMI as objective functions. Higher proton statistics generate comparable results to the lower proton statistics. In particular, in integration-mode mean the higher proton statistics show slightly faster convergence of the NMI than the lower ones. However, this does not translate into better accuracy.
Integration-mode max has demonstrated better performance than integration-mode mean, but remains inferior to list-mode. List-mode demonstrates faster convergence, but also instability of the objective function with the consideration of two proton radiographies. The stability of the results is affected by the limited amount of information provided by two proton radiographies. Different to integration-mode, the multiplicity and accuracy of the information contained in list-mode radiographies limit the searching space of the 2D-3D DIR. In general, the anatomical changes (due to breathing) are approximately compensated with two orthogonal proton radiographies and minor improvement is observed with nine proton radiographies. The accuracy of comparable proton statistics, but different number of proton radiographies (uniformly covering 180 • ) for the NMI as the objective function are compared in figure 8. As expected, comparable results are obtained. Figure 7. Overlay of the ground truth pCT and the calibrated treatment planning CT (anatomical changes in pink/ green), prior to and after 2D-3D DIR and comparing the RMSE and the NMI as objective functions in anthropomorphic phantom for the analytical simulation, nine list-mode radiographies (A), nine integration-mode max radiographies (B) and nine integration-mode mean radiographies (C). Highlighted region corresponds to 100% proton statistics (in red) and 25% proton statistics (in orange).

Monte Carlo simulation of clinical data
The comparison between list-mode, integration-mode max and mean is reported in figures 9 and 10, for different number of proton radiographies. The second stage for integration-mode max and mean (figure 10) is indicated to enable fair comparison with list-mode (figure 9). The selection of two orthogonal list-mode proton radiographies generates comparable, but better results than the selection of nine proton radiographies from integration-mode max. The superiority of integration-mode max over mean is observed. The minimum RMSE of RSP achieved at two stages with two list-mode radiographies is not achieved at three stages with nine integration-mode radiographies. In figure 11, the overlay of the ground truth pCT and the calibrated treatment planning CT, prior to and after 2D-3D DIR, is shown for both list-mode and integration-mode proton radiographies.

Discussion
In geometrical simulation, the achievable accuracy with a minimum of nine proton radiographies uniformly covering the 180 • (RMSE of RSP equal to 0.0593) results adequately for promising ART in proton therapy, being comparable to conventional 3D-3D DIR (RMSE of RSP equal to 0.0536). However, the achieved accuracy remains inferior to the one obtained with conventional 3D-3D DIR, also for a large number (180) of proton radiographies (RMSE of RSP equal to 0.0585). The computational time plays a major role in the exploration of registration parameters. Thus, the finer grid spacing is limited to the geometrical simulation. Therefore, random sampling for objective function calculation is considered for computation time optimization (Klein et al 2007). In 3D-3D DIR, the choice of NMI as the objective function is prevented by the RSP discretization of the anthropomorphic phantom that compromises the numerical convergence of the optimization algorithm. According to the chosen metrics in image domain (RMSE of RSP and Dice similarity coefficients), the results suggest the need for stopping criteria. A deterioration of the metric is noted as the number of iterations increases, especially if the RMSE is chosen as the objective function (figures 4, 5, 6 and 10). However, the stopping criteria have to be introduced, relying on a critical evaluation of the adopted metrics, especially when contradictory deteriorations of different metrics are observed (e.g. increment of Dice similarity coefficients while reduction of the RMSE of RSP, or vice versa). Therefore, the development of stopping criteria in combination with different metrics in image domain is ongoing. In particular, the consideration of the objective function in projection domain, typically stable when the metric in image domain starts deteriorating, can be combined with the evaluation of the Jacobian of the deformation field. A preliminary assessment demonstrates that features (size and connectivity) of the negative regions of the Jacobian are correlated to the instability of the results. Eventually, the evaluation of the accuracy of the treatment planning based on the calibrated treatment planning CT, compensated for the anatomical changes, is anticipated.
The accuracy of the 2D-3D DIR strongly depends on the applied anatomical changes. In this regard, the adopted anthropomorphic phantom represents the worst-case scenario (i.e. large and complex anatomical changes extending across the entire image). Despite the large and complex anatomical changes, the two major components of the deformation field lay on the lateral (i.e. chest wall expansion) and frontal (i.e. diaphragm motion) radiographies. Depending on the amount of information provided by the radiographies as well as the accuracy of this information, a minimum number of radiographies is defined to achieve the Figure 11. Overlay of the ground truth pCT and the calibrated treatment planning CT, prior to and after 2D-3D DIR and comparing the RMSE and the NMI as objective functions in clinical data, for two list-mode radiographies (A), nine integration-mode max radiographies (B) and two integration-mode max radiographies (C). Highlighted in the first column is the Gaussian component of the deformation field (top), deformation field overlaid (in green) to the ground truth pCT (center) and applied anatomical changes (in pink/ green) overlaid to the ground truth pCT (bottom). optimal compensation of anatomical changes. In geometrical simulation, the amount of information is not redundant and is highly accurate. Nine radiographies are sufficient to optimally capture the anatomical changes. When the same amount of information is affected by inaccuracies as in integration-mode radiographies, the achievement of the optimal compensation is hindered, so that nine radiographies perform similar to two orthogonal radiographies. In the presence of a redundant amount of information as in list-mode radiographies, the comparable performance between two and nine radiographies comes with the disadvantage of the instability of the objective function for two radiographies. For list-mode radiographies, the inaccuracies due to the anatomical changes in the proton DRRs (based on the estimated proton trajectories) are dominant when compared to the inaccuracies related to the estimated proton trajectory under the assumption of uniform water. In analytical simulation of the anthropomorphic phantom, proton statistics play a minor role.
Monte Carlo simulation of proton radiographies based on clinical data proved that reasonable proton statistics for 2D-3D DIR correspond to clinically acceptable imaging dose . In agreement with previous studies (Schulte et al 2005, Murphy et al 2007, the imaging dose of realistic clinical proton tomography obtained with the results of 180 radiographies is much lower (~2 mGy) than the one of commercial x-ray tomography . The artificial deformation field applied to clinical data suggests the potentiality of compensating localized anatomical changes at clinically acceptable imaging dose. Generally, list-mode performs better than integration-mode, but the stability of the results depends on the number of proton radiographies that determine the amount of information about the anatomical changes. Given the same number of proton radiographies, integration-mode benefits from the broadened searching space of the 2D-3D DIR provided by the reduced amount of information and reduced accuracy of such information, especially at anatomical interfaces. Thus, integration-mode max performs better than integration-mode mean. Optimization of integration-mode is proposed when the entire WET histogram for each pencil beam is exploited (Deffet et al 2017;Seller Oria et al 2018). Therefore, the forward-projection of the calibrated treatment planning CT requires dedicated calculation (simulation), different to the straight central pencil beam axis.
Future studies will include 2D-3D DIR in the presence of realistic anatomical changes combined with rigid patient position inaccuracies. In particular, in-room CBCT images acquired prior to treatment delivery will be used to estimate the corresponding realistic (instead of artificial) deformation field. The anatomical changes are expected to be accurately compensated, relying on proton radiographies. Information about patient position and anatomical changes will be retrieved by proton radiographies, similar to x-ray radiographic and tomographic imaging. In addition, information about RSP will be also derived, thus enabling the combination of the 2D-3D DIR with the optimization of the empirical calibration of the treatment planning CT directly based on the adopted proton radiographies (Schneider et al 2005). On the one hand, the optimization of the calibration will require adequate registration between proton radiography and treatment planning CT. On the other hand, the registration will benefit from the forward-projection of an accurately calibrated treatment planning CT. The framework will therefore consider either sequential or joint implementation of the two methodologies. In this investigation, the model of realistic detectors will be embedded in Monte Carlo simulation. Eventually, the use of a limited number of proton radiographies for both adaptation and calibration of the treatment planning CT is embedded within the same framework, thus making the most of the native in-room radiographic imaging in proton therapy.

Conclusion
This work first investigates a 2D-3D DIR using a limited number of proton radiographies to compensate the anatomical changes in the treatment delivery scenario with respect to the treatment planning CT. Anatomical changes in an anthropomorphic phantom and clinical data are considered, relying on geometrical and analytical simulations as well as Monte Carlo simulation, to differently control physical and geometrical parameters. List-mode and integration-mode proton radiographies with different proton statistics are simulated to compare the performance of different detector configurations. The results show that the 2D-3D DIR is capable of compensating the considered anatomical changes. The minimum number of radiographies depends on the amount and complexity of the anatomical changes and the adopted detector configurations, but less on proton statistics, thus allowing the usage of low dose proton radiographies in 2D-3D DIR. Methodological improvements as well as computational optimization are necessary towards the clinical implementation of the proposed 2D-3D DIR. Eventually, the accuracy of the treatment (re)planning based on the calibrated treatment planning CT, compensated for the anatomical changes, will be investigated for ART in proton therapy based on a limited number of proton radiographies. Figure A3. Dice similarity coefficients for segmented anatomical region of interest heart (A), lesion (B), liver (C) and lung (D) prior to and after 2D-3D DIR and 3D-3D DIR for different proton statistics ('pS') on analytical simulation of integration-mode max (stages 1 and 2 correspond to 64 and 32 grid spacing, respectively). Figure A4. Dice similarity coefficients for segmented anatomical region of interest heart (A), lesion (B), liver (C) and lung (D) prior to and after 2D-3D DIR and 3D-3D DIR for different proton statistics ('pS') on analytical simulation of list-mode (stages 1 and 2 correspond to 64 and 32 grid spacing, respectively).