Radiation therapy is a local treatment option for small hepatocellular carcinoma (HCC) tumors in inoperable patients with a good liver function (Child-Pugh A) (1). However, HCC patients often present with a considerable tumor burden and an underlying cirrhotic liver and even low doses of radiation to the liver leads to a high risk of developing radiation induced liver disease (RILD) in patients with cirrhosis (2). Since RILD is a severe condition that can lead to liver failure and death, it is crucial to minimize the dose to the normal liver tissue surrounding the tumor (3).

Compared to photon based radiation therapy, pencil beam scanning (PBS) proton therapy can often provide more conformal target dose distributions with less healthy tissue irradiation (4, 5). Proton therapy is therefore increasingly used in the treatment of HCC (6, 7). However, liver tumors often exhibit large and variable respiratory motion during treatment (8), which can cause considerable deviations between the delivered and planned doses. Due to interplay effects and a high sensitivity to water equivalent path length changes, PBS proton therapy is particularly prone to dosimetric uncertainties caused by target motion (9–11) and international guidelines underline the special need for motion management in PBS proton therapy (12, 13). Hence, respiratory gating, where the beam is only turned on during specific phases of the breathing cycle has been proposed and implemented in proton therapy to mitigate tumor motion effects (14–17). Still, residual motion within the gating window is of concern (16, 18).

Reconstruction of the actual delivered tumor dose at a fraction requires knowledge of the internal motion during treatment delivery and synchronization of this motion with the beam delivery. One method is to calculate the dynamic 4D dose (D4DD) by ascribing a specific phase of a 4DCT scan to each delivered spot, use this to calculate phase-specific doses in each 4DCT phase and accumulate these doses in a reference phase by deformable image registration (12, 13). This method has been implemented clinically for PBS carbon therapy (19) and proton therapy (20, 21) using a waist belt for respiratory monitoring during beam delivery. A limitation of the D4DD is that it neglects setup errors and assumes that the internal motion during treatment is well described by the respiration signal and identical to the motion in the 4DCT. However, liver motion is known to be highly variable and often poorly represented by 4DCT scans that by nature only capture one (random) respiratory cycle at each longitudinal position within the patient (22, 23). To overcome these limitations, Yamada et al. monitored the internal motion of implanted fiducial markers in the liver during gated proton PBS delivery by a gantry-mounted stereoscopic fluoroscopic x-ray imager (24). By combining the internal motion with the spot delivery timing in machine log files the authors reconstructed the tumor dose by a spot shift method that can account for arbitrary rigid motion (25). However, although many modern conventional proton facilities are equipped with dual x-ray imagers, these can typically not be used during treatment delivery. Consequently, target motion monitoring during treatment is normally not available even though it is recommended by international guidelines for proton therapy of moving targets (12, 13).

In this study, we introduce a method to overcome the limitations of conventional proton facilities in internal tumor motion monitoring during proton PBS treatment. The method uses an external-internal correlation model (ECM) to estimate the internal tumor motion from an external respiratory signal and combines the internal motion with spot delivery timing in machine log files to estimate the tumor position during each spot delivery. The motion is then combined with the spot shift dose reconstruction method to estimate the tumor doses for individual HCC treatment fractions.

Figure 3 shows an example of the tumor motion trajectory during the setup CBCT and the post-treatment CBCT at one fraction for patient 3. At this fraction, the online registration errors were 0.4-0.8 mm ( Figure 3 ). The mean online registration error was in general sub-millimeter for all patients ( Table 1 ). The example case in Figure 3 had a small cranial and posterior drift of the tumor exhale position from setup CBCT to post-treatment CBCT (black arrows). Averaged over all fractions and patients a similar trend was seen with mean ± standard deviation (SD) drift motion of 0.0 mm ± 0.8 mm (LR), 1.3 mm ± 1.3 mm (CC), and -0.7 mm ± 1.0 mm (AP) ( Table 1 ).

The example ECM presented in Figure 3 (blue curves) had an accuracy close to the mean RMS fit error for patient 3, while patients 1 and 2 had larger RMS fit errors up to 2.1 mm ( Table 1 ). Over all patients, the mean ( ± SD) RMS fit error of the ECM was 0.5 mm ± 0.4 mm (LR), 1.5 mm ± 0.8 mm (CC), and 1.0 mm ± 0.6 mm (AP).

Figure 4 presents a typical example of the RPM signal and the ECM estimated tumor motion at a fraction, synchronized with the spot delivery times. Due to the gating latency the beam started 240 ms into the gating window and continued 80 ms after the RPM block moved out of the gating window ( Figure 4C ). The treatment error is reported in Table 1 for each patient. Over all delivered spots the RMS treatment error was 1.3 mm (LR), 2.8 mm (CC), and 1.7 mm (AP), while the mean ( ± SD) 3D treatment error per patient was 3.9 mm ± 1.9 mm (patient 1), 3.7 mm ± 0.6 mm (patient 2) and 2.6 mm ± 1.7 mm (patient 3).

The maximum tumor motion range during a fraction was 6.4 mm (LR), 27.9 mm (CC), and 19.2 mm (AP) during field delivery independent of beam-on status and 2.7 mm (LR), 10 mm (CC), and 7.1 mm (AP) during beam-on periods. The mean tumor motion range during a single fraction was usually more than halved with gating compared to the full motion range ( Table 1 ).

Large dose deterioration occurred at single fractions due to interplay effects with D2% being up to 4.7%-points higher than planned and D98% up to 4.4%-points lower than planned ( Figures 5 , 6 ; Table 2 ). After 5-6 fractions the interplay effects tended to smear out due to averaging effects such that D2% and D98% converged towards the planned values ( Figure 6 ). On average the CTV HI was larger than planned by 3.7%-points (range: 1.0-6.6%-points) for individual fractions and by 0.7%-points (range: 0.3-1.1%-points) for the average dose of 5 or 6 fractions ( Table 2 ). The mean dose to the healthy liver tissue, averaged over all analyzed fractions, was different from the planned dose by 0.3%-points (patient 1), 1%-points (patient 2) and -0.1%-points (patient 3).

With the present study we have developed and clinically demonstrated a method to estimate the internal target motion and its consequence on dose delivery during proton therapy of liver cancer at a standard equipped proton facility. As pointed out by recent international guidelines such monitoring is important for PBS proton therapy of moving targets but typically not commercially available (12, 13). The motion estimation was based on an ECM that was constructed at each treatment fraction using external RPM motion synchronized with internal 3D tumor motion extracted from CBCT projections. An in-house developed method that has been validated on a group of ten patients (11, 25) was subsequently applied to reconstruct the motion-including CTV dose. Considerable interplay effects at single fractions tended to smear out after more fractions.

Additionally, we investigated the accuracy of the online CBCT match to determine the exhale marker positions. Due to motion smearing and motion artifacts, the manual online match is subjective and prone to human errors, while the estimated marker trajectories provided an objective measure of the exhale position during CBCT. However, with differences between the online manual registration and offline trajectory-based marker match close to the resolution of the CT scan (2mm CC, ~1mm in-plane), the online match accuracy was acceptable. Yet, at a few individual fractions, slightly larger discrepancies in the CC direction were observed (up to 3.2 mm). Offline inspection of the online match revealed that these discrepancies could be ascribed to the online procedure being prone to human subjectivity and performed under time pressure with the patient waiting for treatment.

During CBCT scan and treatment delivery, large motion variations between individual respiratory cycles and total motion amplitudes of 2-3 cm, most prominent in the CC direction, were observed ( Figures 3 , 4 ). Such motion is on par with previous studies of internal motion during radiotherapy of tumors in the liver (8, 32–35). Due to the extended time typically spent near the exhale phase of the respiratory cycle, limiting the treatment to an approximate 50% duty cycle around exhale more than halved the motion during beam-on and also reduced motion variation between treatment fractions ( Table 1 ). The resulting mean geometrical treatment errors during proton spot delivery were thereby also limited to a few millimeters, including errors introduced by the rather small baseline shift of 0 - 3 mm typically observed between setup CBCTs and post-treatment CBCTs ( Table 1 ). Notably, though not observed in the limited cohort of this study, baseline shifts and resulting treatment errors of higher magnitude can be expected during liver treatments of some patients (8, 33). Still, the initial findings of the present study confirmed the usability of externally guided respiratory gating, which has also been proposed by other groups for reducing the internal motion during liver proton therapy delivery (14–16).

Despite the use of respiraory gating, our dose reconstructions showed considereable dose deteriorations on a single-fraction level caused by interplay between proton spot delivery and target motion ( Figures 5 , 6 ). This finding clinically confirms the simulation results of Zhang et al. who concluded that respiratory gating alone was insufficient to mitigate interplay effects in PBS proton therapy (16). However, fractionation tends to reduce interplay effects by averaging out local over and under dosage over several treatment fractions (36–39). In the present study, the interplay effects were almost averaged out after 5-6 fractions ( Figure 5 ). Nevertheless, for hypofractionation treatments, one may need to combine gating with repainting (16, 40).

A few previous clinical studies investigated the dosimetric consequences of respiratory motion during particle therapy. Richter et al. (19) and Meijers et al. (21) combined machine log files with the spot delivery with an external respiratory signal obtained during treatment and used this to calculate the D4DD by distributing each spot delivery into corresponding phases of 4DCT scans. These studies also found that fractionation effectively mitigated interplay effects. For tumors in the thoraic region, anatomical changes such as presence of fluid or tumor shrinkage, caused more severe dosimetric changes (21). Limitations of the D4DD include neglectance of setup errors, the assumption of identical anatomy and respiratory motion amplitude at treatment as in the 4DCT and the dependence on deformable image registration for dose accumulation in a reference 4DCT phase. An advantage of the 4DCT based dose reconstruction, compared to our method, is that it includes estimations of dose degration caused by internal anatomical changes between individual 4DCT phases. Since our spot shift dose reconstruction was based on a single phase of the planning 4DCT it only considers the effects of rigid intrafraction motion. A potential improvement could be a hybrid dose reconstruction method that extends phase specific dose calculations in each 4DCT phase with spot shifts that accounts for the tumor motion during treatment that goes beyond the motion in the 4DCT scan. Such an extension of our method to individual 4DCT phases would improve the D4DD reconstruction method recommended by international guidelines to also include actual intrafraction motion and not only motion observed in 4DCTs (12, 13). However, a persisting challenge with such 4DCT dose accumulation is the reliance on deformable image registration and dose warping between CT scans, which is a procedure with considerable uncertainty (41).

In a recent and closely related study by Yamada et al. at Hokkaido University, machine log files were combined with intra-treatment monitoring of internal fiducial markers by stereoscopic x-ray fluoroscopy during liver proton therapy (24). The same dose reconstruction method as applied in the present study was used (25). Respiratory exhale gating guided by direct internal tumor motion monitoring was a big advantage of the study by Yamada et al. compared to ours. Tight internal gating windows of ±2 mm along each direction may also explain why Yamada et al. found considerably smaller mean 3D tumor position errors during spot delivery (0.8-1.3 mm) than the present study (2.6-3.9 mm) where only external monitoring with a gating window corresponding to approximately 50% duty cycle was applied. For comparison, a previous study on liver SBRT with internal exhale gating based on implanted electromagnetic markers and gating windows of ±3 mm (LR/AP) and ±4 mm (CC) found mean 3D tumor position errors of 1.2-3.0 mm (8). Despite the larger treatment errors in the present study the estimated delivered CTV dose was close to the planned dose when averaged over 5-6 fractions. Hence, the robust SFUD planning approach combined with CBCT-based setup to internal markers and fractionated exhale gated treatments provided appropriate mitigation of intrafraction motion. Further studies including more patients and the effects of inter-fraction deformations (e.g., by dose reconstruction on weekly control 4DCT scans) are necessary to conclude on the overall treatment quality. For example, deformations may affect the proton range and the assumption that implanted markers serve as an accurate surrogate for the CTV position, especially if they are not implanted near the tumor (42). Ultimately, comprehensive all-inclusive fraction-specific dose reconstruction and dose accumulation could be used to trigger plan-adaptations in case of unacceptable dose coverage.

A limitation of the current study is the indirect estimation of the tumor motion during treatment by an ECM. Although an ECM of the day was built, the accuracy cannot be expected to be better than the ECM fit, which had mean RMS errors of almost half of the treatment errors and exceeding 2 mm in the CC direction for one patient ( Table 1 ). Furthermore, high intrafraction stability of the ECM with the ability to detect internal baseline shift cannot in general be assumed (43–46). Optimally, the ECM should be built based on the setup CBCT and then validated by the post-treatment CBCT. However, phantom tests showed discrepancies between couch shifts and changes in the RPM signal that hindered correct adjustment of an ECM from a setup CBCT to allow usage after setup couch corrections. For this reason, intra-treatment motion could only be investigated for fractions with a post-treatment CBCT in this study. The ECM stability issue may be addressed by capturing a series of stereoscopic images before each field delivery to confirm and update the ECM.

It is worth noting that the RPM system available for proton therapy clearly lags behind the optical monitoring system available from the same vendor for photon radiotherapy (TrueBeam, Varian) with a stereoscopic camera that allows adaptation of the ECM to couch shifts (47). Furthermore, the cumbersome manual synchronizations of the RPM log files with CBCT projections and treatment machine log files in the current study ( Figure 1 ) would not be needed if the respiratory monitoring, imaging, and beam delivery systems were similarly well-integrated as on a TrueBeam linear accelerator. The manual synchronization of RPM log files with machine log files was only possible in this study because the beam pauses in the machine log files had a unique temporal pattern that could be matched with the RPM gating signal ( Figure 1E ). This synchronization would not be possible for non-gated treatments. Since the synchronization can only be performed post-treatment, it is currently a barrier for online real-time dose reconstruction, which has been demonstrated clinically for liver SBRT (48), but could be even more relevant for proton PBS.

In summary, dose reconstruction including the effects of setup errors, rigid motion, interplay effects and the smearing hereof after 5-6 fractions was performed for HCC proton therapy. For the included patients, it showed that our treatment strategy of exhale gating resulted in an acceptable CTV dose coverage. With a smoother workflow and automation this could be used to trigger a plan adaptation if the CTV dose coverage turned out to be unacceptable. Since CTV dose deficits could also be caused by interfractional changes, the motion-including dose reconstruction should ideally be extended to account for such changes, for example by applying it on the anatomy of weekly 4DCTs.

A method to estimate internal tumor motion and reconstruct the motion-including fraction dose for PBS proton therapy in the liver was developed and successfully demonstrated clinically at a conventional proton facility.

The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding author.

The studies involving human participants were reviewed and approved by the regional research ethics committee and the regional register of research projects (ref. number 742503, case number 1-16-02-320-21). The research was conducted in accordance with the principles of the Helsinki Declaration and local statutory requirements. The patients/participants provided their written informed consent to participate in this study.

EW, HM, BW, MH, JT, LS, PP designed the treatment and imaging protocol for the study. PP, SN, EW conceived, designed and drafted the manuscript and all authors revised the manuscript. JT, LS, SN, EW collected the data. Software was developed by JB (marker segmentation) and PP (segmentation correction, 3D trajectory estimation, data synchronization, dose reconstruction), and streamlined by SN who performed the data analysis. All authors contributed to the article and approved the submitted version.