Intrafraction tumor motion monitoring and dose reconstruction for liver pencil beam scanning proton therapy

Background Pencil beam scanning (PBS) proton therapy can provide highly conformal target dose distributions and healthy tissue sparing. However, proton therapy of hepatocellular carcinoma (HCC) is prone to dosimetrical uncertainties induced by respiratory motion. This study aims to develop intra-treatment tumor motion monitoring during respiratory gated proton therapy and combine it with motion-including dose reconstruction to estimate the delivered tumor doses for individual HCC treatment fractions. Methods Three HCC-patients were planned to receive 58 GyRBE (n=2) or 67.5 GyRBE (n=1) of exhale respiratory gated PBS proton therapy in 15 fractions. The treatment planning was based on the exhale phase of a 4-dimensional CT scan. Daily setup was based on cone-beam CT (CBCT) imaging of three implanted fiducial markers. An external marker block (RPM) on the patient’s abdomen was used for exhale gating in free breathing. This study was based on 5 fractions (patient 1), 1 fraction (patient 2) and 6 fractions (patient 3) where a post-treatment control CBCT was available. After treatment, segmented 2D marker positions in the post-treatment CBCT projections provided the estimated 3D motion trajectory during the CBCT by a probability-based method. An external-internal correlation model (ECM) that estimated the tumor motion from the RPM motion was built from the synchronized RPM signal and marker motion in the CBCT. The ECM was then used to estimate intra-treatment tumor motion. Finally, the motion-including CTV dose was estimated using a dose reconstruction method that emulates tumor motion in beam’s eye view as lateral spot shifts and in-depth motion as changes in the proton beam energy. The CTV homogeneity index (HI) The CTV homogeneity index (HI) was calculated as D2% − D98%D50% ×100% . Results The tumor position during spot delivery had a root-mean-square error of 1.3 mm in left-right, 2.8 mm in cranio-caudal and 1.7 mm in anterior-posterior directions compared to the planned position. 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. Conclusions A method to estimate internal tumor motion and reconstruct the motion-including fraction dose for PBS proton therapy of HCC was developed and demonstrated successfully clinically.

Background: Pencil beam scanning (PBS) proton therapy can provide highly conformal target dose distributions and healthy tissue sparing. However, proton therapy of hepatocellular carcinoma (HCC) is prone to dosimetrical uncertainties induced by respiratory motion. This study aims to develop intra-treatment tumor motion monitoring during respiratory gated proton therapy and combine it with motion-including dose reconstruction to estimate the delivered tumor doses for individual HCC treatment fractions.
Methods: Three HCC-patients were planned to receive 58 GyRBE (n=2) or 67.5 GyRBE (n=1) of exhale respiratory gated PBS proton therapy in 15 fractions. The treatment planning was based on the exhale phase of a 4-dimensional CT scan. Daily setup was based on cone-beam CT (CBCT) imaging of three implanted fiducial markers. An external marker block (RPM) on the patient's abdomen was used for exhale gating in free breathing. This study was based on 5 fractions (patient 1), 1 fraction (patient 2) and 6 fractions (patient 3) where a post-treatment control CBCT was available. After treatment, segmented 2D marker positions in the post-treatment CBCT projections provided the estimated 3D motion trajectory during the CBCT by a probability-based method. An externalinternal correlation model (ECM) that estimated the tumor motion from the RPM motion was built from the synchronized RPM signal and marker motion in the CBCT. The ECM was then used to estimate intra-treatment tumor motion. Finally, the motion-including CTV dose was estimated using a dose reconstruction method that emulates tumor motion in beam's eye view as lateral spot shifts and in-depth motion as changes in the proton beam energy. The CTV homogeneity index (HI) The CTV homogeneity index (HI) was calculated as D2% − D98 % D50 % Â 100 %.
Results: The tumor position during spot delivery had a root-mean-square error of 1.3 mm in left-right, 2.8 mm in cranio-caudal and 1.7 mm in anterior-posterior directions compared to the planned position. On average, the CTV HI was larger

Introduction
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)(15)(16)(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.

Patients and treatment planning
Three patients with HCC underwent proton PBS in April-September 2022 in a national phase II clinical trial that allowed inclusion of both large tumors and Child-Pugh B patients. The trial was approved by the relevant ethics committees (ClinicalTrials.gov Identifier: NCT05203120). Three gold or platinum fiducial markers with dimensions of 0.75 mm × 5 mm (Visicoil ™ ) were implanted near the tumor the day before planning CT scanning. An internal clinical target volume (iCTV) was formed as the union of the CTV in the five phases of a 10-phase 4-dimensional CT scan (4DCT) that were closest to full exhale. It corresponded to exhale respiratory gating with approximately 50% duty cycle. A 3-field proton plan was created on the exhale phase of the 4DCT using a commercial treatment planning system (TPS, Eclipse 16.01.10, Varian, a Siemens Healthineers Company, Palo Alto, CA) and dose calculation algorithm (Varian Proton Convolution Superposition 16.1.0). Robust single field uniform dose (SFUD) optimization was performed with ±4.5% range uncertainty and ±5 mm shifts in the left-right (LR) and anterior-posterior (AP) directions and ±7 mm shifts in the cranio-caudal (CC) direction. The treatment plans used beam energies in the range 71-153 MeV, for which the spot size in air is 4-6mm (1 standard deviation). Each field had 528-2134 spots and a total of 6617-16047 monitor units. The prescribed mean iCTV dose was 58 GyRBE (Patient 1 and 3) for central tumors (≤2 cm from porta hepatis) or 67.5 GyRBE (Patient 2) for peripheral tumors (>2 cm from porta hepatis) in 15 fractions.

Treatment delivery and imaging
Daily patient setup was based on a free-breathing CBCT scan in which the estimated exhale positions of the motion-blurred fiducial markers were matched with the planning CT. CBCT scan was done using Standard ProBeam CBCT imaging system with Paxscan 4030CB flat panel detectors. The resolution of the image detector was 0.388 mm/pixel in both directions with source-to-imager distance (SID) of 3700 mm and source-to-axis distance (SAD) of 2700 mm. During the CBCT acquisition and throughout the whole treatment session the position of a marker block (Real-time Position Management System, RPM, Varian) on the patient's abdomen was recorded with a camera. During treatment the RPM signal was used for respiratory gating with a gating window adjusted before treatment to obtain a duty cycle of approximately 50% centered around the exhale phase in accordance with the iCTV construction. A post-treatment control CBCT scan was captured at 6, 1, and 7 fractions for patients 1, 2, and 3, respectively. The RPM log file was missing for one of these treatment fractions for patient 1 and patient 3. The analysis presented in this study requires a post-treatment CBCT and an RPM log file and was therefore only made for 5, 1 and 6 fractions for patients 1, 2, and 3, respectively.

Data analysis
After the treatments the fiducial markers were segmented in each raw 2D CBCT-projection (~1000 images per CBCT) using an automatic method (26) followed by manual inspection and semiautomatic correction of failed segmentations. The 3D motion trajectory of each marker during CBCT was estimated by a probability-based method (27) and the marker group centroid motion was used as a surrogate for the tumor motion. The exhale period was defined as the time within the 95 th -100 th percentile of the markers position in the CC direction for each CBCT. This was used to determine the exhale position in each direction of motion as the mean marker position during the exhale period. For the setup CBCT scans, the resulting exhale tumor position was used to determine the optimal setup couch correction for marker alignment with the planned marker positions in exhale. This is similar to the trajectory-based setup introduced for non-gated treatments in (28). The online registration error was then calculated as the difference between the retrospective trajectory-based patient setup and the actual couch correction based on online 3D/3D registration of the setup CBCT with the planning CT scan. Furthermore, the intrafraction baseline drift of the exhale position between the setup CBCT and the post-treatment CBCT was determined as the difference between their respective exhale positions.
The analysis in this study required synchronization of the RPM signal with the projection images of the post-treatment CBCT (to establish an ECM) and with the delivery time of each proton spot (to perform dose reconstruction). Synchronization between RPM and CBCT projections was obtained by placing a 3 mm diameter tungsten sphere on the RPM block such that it was visible in most of the CBCT projections ( Figures 1A, B). After treatment the 3D motion of the tungsten sphere during the CBCT scan was estimated from its projected motion in the CBCT projections (27), and its AP motion was temporally aligned with the RPM motion in gating log files ( Figure 1C). This synchronization provided the logged RPM position at the acquisition time of each CBCT projection.
The synchronization between RPM and spot delivery times was based on the gate-on signal in the gating log files. The logged gate-on signal specifies the time intervals when the RPM block is inside the gating window, but it does not account for the gate-on latency between entering the gating window and beam-on and the gate-off latency between exiting the gating window and beam-off. The gate-open times during which the beam could potentially be turned on were estimated from the logged gate-on signal by assuming a gate-on latency of 240 ms and a gate-off latency of 80 ms ( Figure 1D). These latencies were measured using the method proposed by Worm et al. (29) and rounded to an integer number of gating log file samples (40ms resolution). Next, a comparison of the gate-open times with the actual spot delivery times in machine log files (30) provided the synchronization between RPM log files and spot delivery times ( Figure 1E). While the machine log files specified the duration of each spot with microsecond resolution it did not directly specify the beam-off times occurring during larger spots shifts, energy shifts and gate-off periods. However, by using the logged number of User Datagram Protocol (UDP) messages received between each spot the beam-off times were estimated with a scaling factor uncertainty of a few percent. During the synchronization of the machine log files with the gating log files the beam-off times were scaled to fit the time scale in the gating log file.
To estimate the tumor motion during treatment delivery an augmented linear ECM that estimated the tumor motion during the post-treatment CBCT from the synchronized RPM motion (31) was built (Label 1 in Figure 2): Here, INT and EXT are internal tumor and external RPM motion as a function of time (t). The coefficients A, B and C and the time delay t are fitting parameters. The augmentation term B.EXT (t-t) accounts for hysteresis and phase differences between internal and external motion (31). A, B and C were optimized individually for each motion direction with least-square fitting while the same value of t was used for all three motion directions.
Next, The ECM was used to estimate the tumor motion throughout the treatment delivery from the intra-treatment RPM signal (Label 2). This synchronization resulted in the ECM estimated tumor position at the time of each spot delivery (Label 3). For each spot, the geometrical treatment accuracy was determined as the tumor position relative to the planned position. The tumor motion range was calculated for each fraction as the difference between the 98 th and the 2 nd percentiles of the tumor position over all three fields during beam-on periods and, for comparison, over the full field durations including both beam-on and beam-off periods. Finally, the CTV dose and the dose to the healthy liver tissue of each analyzed fraction was estimated by a motion-including dose reconstruction method that emulates tumor motion in beam's eye view by shifting each spot in the opposite direction of the tumor motion and in-depth tumor motion as changes in the proton beam energy (25). A motion-encoded plan with these manipulations of the original spot positions and energies was created by an in-house developed Matlab program (Label 4) and imported into the TPS for recalculation (Label 5). The reconstructed CTV doses of each fraction and averaged over all analyzed fractions were compared with the planned doses using the metrics of CTV D98% and D2% (minimum dose received by 98% and 2% of the CTV volume) and the homogeneity index HI% = D2% − D98 % D50 % Â 100 %. Furthermore, the mean healthy liver tissue doses averaged over all analyzed fractions were compared with the planned doses. 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).

Results
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).

Discussion
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)(33)(34)(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 posttreatment 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)(15)(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)(37)(38)(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 planadaptations 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)(44)(45)(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 wellintegrated 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 realtime 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.

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

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

Ethics statement
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.

Author contributions
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.

Conflict of interest
PP is co-inventor on a patent on the 3D motion estimation method applied in this study.
The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.