Continuous time-resolved estimated synthetic 4D-CTs for dose reconstruction of lung tumor treatments at a 0.35 T MR-linac

Objective. To experimentally validate a method to create continuous time-resolved estimated synthetic 4D-computed tomography datasets (tresCTs) based on orthogonal cine MRI data for lung cancer treatments at a magnetic resonance imaging (MRI) guided linear accelerator (MR-linac). Approach. A breathing porcine lung phantom was scanned at a CT scanner and 0.35 T MR-linac. Orthogonal cine MRI series (sagittal/coronal orientation) at 7.3 Hz, intersecting tumor-mimicking gelatin nodules, were deformably registered to mid-exhale 3D-CT and 3D-MRI datasets. The time-resolved deformation vector fields were extrapolated to 3D and applied to a reference synthetic 3D-CT image (sCTref), while accounting for breathing phase-dependent lung density variations, to create 82 s long tresCTs at 3.65 Hz. Ten tresCTs were created for ten tracked nodules with different motion patterns in two lungs. For each dataset, a treatment plan was created on the mid-exhale phase of a measured ground truth (GT) respiratory-correlated 4D-CT dataset with the tracked nodule as gross tumor volume (GTV). Each plan was recalculated on the GT 4D-CT, randomly sampled tresCT, and static sCTref images. Dose distributions for corresponding breathing phases were compared in gamma (2%/2 mm) and dose–volume histogram (DVH) parameter analyses. Main results. The mean gamma pass rate between all tresCT and GT 4D-CT dose distributions was 98.6%. The mean absolute relative deviations of the tresCT with respect to GT DVH parameters were 1.9%, 1.0%, and 1.4% for the GTV D 98%, D 50%, and D 2%, respectively, 1.0% for the remaining nodules D 50%, and 1.5% for the lung V 20Gy. The gamma pass rate for the tresCTs was significantly larger (p < 0.01), and the GTV D 50% deviations with respect to the GT were significantly smaller (p < 0.01) than for the sCTref. Significance. The results suggest that tresCTs could be valuable for time-resolved reconstruction and intrafractional accumulation of the dose to the GTV for lung cancer patients treated at MR-linacs in the future.


Introduction
Stereotactic magnetic resonance imaging-guided radiotherapy (MRgRT) of peripheral and central lung tumors has been clinically implemented at MR-guided linear accelerators (MR-linacs) over the last decade.First clinical reports show that high ablative doses can safely be delivered to lung tumors with low rates of high-grade toxicity, even to challenging entities such as central or ultracentral lung tumors (Finazzi et al 2019, Henke et al 2019, Regnery et al 2022, 2023) in few (3-12) fractions (Finazzi et al 2020a, 2020b, Crockett et al 2021) or even just a single fraction (Finazzi et al 2020c, Palacios et al 2022).
Online treatment plan adaptation before each treatment fraction based on the daily patient anatomy can improve target coverage and organ-at-risk (OAR) sparing for MRgRT of lung cancer (Finazzi et al 2019, Nierer et al 2022).In clinical practice at MR-linacs today, the goal of online plan adaptation is to restore the original treatment plan quality in terms of planned target coverage and OAR sparing in presence of interfractional changes, without explicit consideration of the accumulated dose that has been delivered in the previous treatment fractions (Finazzi et al 2019, Klüter 2019).Interfractional dose accumulation based on fractional magnetic resonance imaging (MRI) datasets and daily delivered dose distributions at MR-linacs is being investigated for a more accurate assessment of the actually delivered dose in research studies for different treatment sites, including the liver (Wahlstedt et al 2023) and lung (Rabe et al 2023).
For an even more accurate reconstruction of the delivered dose, the effects of intrafractional changes including respiratory-induced tumor and OAR motion must be considered.This is of particular importance for ultra-hypofractionated treatment regimes (Finazzi et al 2020b) or single fraction treatments (Finazzi et al 2020c, Palacios et al 2022) for which dose inhomogeneities do not average out over the course of treatment.Recently, Xiong et al (2022) and Wahlstedt et al (2022) derived first-order approximations of the intrafractionally accumulated dose delivered in online adaptive MRgRT to the prostate and close-by OARs in presence of residual target motion within a gating window.They neglected interplay effects and assumed the shift-invariance of dose distributions (McCarter andBeckham 2000, Sharma et al 2012) in photon therapy and convolved a static dose distribution with the target trajectories observed on planar sagittal cine MR images acquired at MRIdian MRlinacs (ViewRay Inc., Oakwood Village, Ohio, USA).These approximations might be justified for treatment sites with homogeneous tissue densities and small and slow target motion like in the pelvic region.However, larger discrepancies to the actually delivered dose are to be expected when this dose reconstruction method is applied to heterogeneous treatment sites affected by more complex and irregular motion patterns and potentially larger out-of-plane motion (in left-right direction) like lung tumors.
More accurate intrafractional dose accumulation based on real-time MRI requires time-resolved synthetic 4D-computed tomography (CT) data and linac log files (Paganelli et al 2018a, Johnstone et al 2018, Menten et al 2020), which are both currently not available at clinical MR-linacs.In the future, the treatment plan could be updated after each fraction for the remaining fractions, while considering the inter-and intrafractionally accumulated dose that was already delivered to the target and OARs.Ultimately, with fast time-resolved synthetic 4D-CT generation methods, dose calculation, and optimization algorithms, the partially delivered dose could be reconstructed and accumulated in real-time during the treatment fraction itself, and the treatment plan could be continuously adapted based on this information (Kontaxis et al 2015, Menten et al 2017).Further applications of such time-resolved synthetic 4D-CTs include gating window optimization (Oh et al 2019), 4D and robust treatment plan optimization and analysis (Heath et al 2009, Meschini et al 2022), and the investigation of interplay effects (Rao et al 2012, von Münchow et al 2022).Furthermore, the dose reconstructed over the whole treatment course could serve as input for clinical dose-response modeling in the post-treatment phase (van Herk et al 2018).
Several research groups have developed methods to create time-resolved or respiratory-correlated synthetic 4D-CT datasets for the thoracic or abdominal region based on 3D-CT and time-resolved or respiratorycorrelated 4D-MRI data.The common approach is to first deformably register a 3D-CT image to one 3D-MRI dataset of a 4D-MRI dataset to create a synthetic 3D-CT image.Then, the deformation vector fields (DVFs) of deformable image registrations (DIRs) between the different breathing phases or time steps of the 4D-MRI dataset are applied to the deformed synthetic 3D-CT image to output a synthetic 4D-CT dataset.However, some methods rely on real-time 4D-MRI sequences (Marx et al 2014, Yang et al 2015, Grimbergen et al 2023) that are not clinically available at MR-linacs and typically provide low spatial resolution.Other methods use a respiratory-correlated 4D-MRI as input data instead (Boye et al 2013, Freedman et al 2019, Paganelli et al 2019, Meschini et al 2020), which is time-consuming to acquire and only represents one average breathing cycle that cannot model the intercyclic breathing variations that occur during lung cancer treatments (McClelland et al 2013).Finally, one method relies on a pre-trained motion model (Müller et al 2018), which also uses a respiratory-correlated 4D-MRI as input data and additionally requires pre-treatment model training.Thus, neither of the proposed methods can be currently applied during beam delivery at MR-linacs.
To address these current limitations, we propose a method to create continuous time-resolved estimated synthetic 4D-CTs (tresCTs) at a clinical MR-linac.The method is based on the propagation method, originally proposed for MRI by Paganelli et al (2018b).It uses orthogonal cine MRI data, with sagittal and coronal slices intersecting the moving tumor, as input data to deform a reference 3D-MRI dataset to create a continuous timeresolved 4D-MRI dataset.The propagation method yielded the highest tumor tracking accuracy in an in silico comparison with similar estimated 4D-MRI generation methods (Paganelli et al 2019).Furthermore, the propagation method to create estimated 4D-MRI datasets has been adapted for use at a MRIdian MR-linac by Rabe et al (2021) and experimentally validated with a porcine lung phantom, which showed that the 3D anatomy in the whole phantom thorax could be estimated with high accuracy.Recently, the method has been extended by Meschini et al (2022) to create time-resolved virtual CTs based on a 3D-CT dataset and orthogonal cine MRI data acquired at a 3.0 T diagnostic scanner for offline treatment robustness evaluation in carbon-ion radiotherapy of pancreatic cancer.
In this study, we investigated the applicability of the propagation method to generate tresCTs permitting dose calculation for lung cancer treatments at a clinical MR-linac.We propose and describe the continuous tresCT generation method for application during MRgRT for dose reconstruction and report the results of experimental validation measurements with a porcine lung phantom at a CT scanner and MRIdian MR-linac.The geometric and dose calculation accuracy of the tresCTs was determined by comparing the tresCT datasets to acquired ground truth (GT) respiratory-correlated 4D-CT datasets exploiting the reproducibility of the breathing phantom.

Experimental setup
The MR-compatible porcine lung phantom artiCHEST (PROdesign GmbH, Heiligkreuzsteinach, Germany) (Biederer and Heller 2003) was used for the experiments in this study.This phantom mimics patient breathing with an ex vivo porcine lung and has been used in a previous experimental validation study to evaluate the geometric accuracy of estimated 4D-MRI datasets created with the propagation method (Rabe et al 2021).The phantom and the lung preparation procedure is described in detail in Rabe et al (2021) and Rabe (2022).The main components of the phantom are two double-walled shells filled with water that surround an ex vivo porcine lung, a water-filled silicone diaphragm located inferiorly to the lung, tubes and hoses connecting the phantom with vacuum and pressure pumps, and a motion control system.Ex vivo porcine lungs can be periodically and reproducibly moved by applying air pressure to a silicone membrane at the inferior end of the diaphragm.
The measurements described below were repeated with different motion patterns for two lungs to create three datasets in total (one dataset for lung 1 and two datasets for lung 2).A heated gelatin-water mixture (0.3 g ml −1 ) was injected into each lung at various locations.Upon contact with the lung tissue, the mixture solidified and formed four (dataset 1, lung 1) or three (datasets 2 and 3, lung 2) structures.These nodules (ten nodules in total for all datasets) served as surrogate target lesions and OARs.Several nodules were injected per lung to test the tresCT method for different positions within the lung.Periodic motion patterns with a breathing frequency of 12 cycles min −1 , varying motion amplitudes, and different baseline pressures (determining the maximum lung inflation state) were set at the control system for the three datasets.
The porcine lung phantom was first scanned with a CT scanner, then transported to, and scanned at an MRlinac (datasets 1 and 2).For dataset 3, the phantom was first scanned at the MR-linac before CT imaging.The breathing motion was paused at the inhale position before transport between the two scanners, and the lung position was marked on the upper phantom shell.The experimental setup was connected to an uninterruptible power supply during transport to constantly evacuate the phantom cavity with the vacuum pumps to ensure a stable lung inflation state.The lung position was checked with the marker positions after transport to identify potential deviations between the lung inflation states at the two scanners.

CT data acquisition
The breathing phantom was imaged with a Toshiba Aquilion LB (Canon Medical Systems, Ōtawara, Japan) CT scanner, used in clinical routine for the acquisition of planning CT images for radiotherapy treatment planning.While the phantom was breathing, the projection data for a respiratory-correlated 4D-CT dataset were acquired.Simultaneous to the image acquisition, a surrogate signal correlated to the breathing phase was recorded with a load cell.The load cell is a component of an abdominal pressure belt respiratory gating system (Anzai Medical Co., Ltd, Tokyo, Japan) and was connected to the pressure hose ventilating the diaphragm with a dedicated adapter.The projections were retrospectively assigned to ten breathing phase bins using the vendor's phasebased sorting algorithm to reconstruct a respiratory-correlated 4D-CT dataset (in-plane resolution: 1.074 × 1.074 mm 2 ; slice thickness: 3 mm; acquisition matrix: 512 × 512; x-ray tube voltage: 120 kV) with ten breathing phases (0%-90% with 10% step size).This respiratory-correlated 4D-CT served as the GT dataset in the geometric and dose calculation analyses described below.

MR-linac data acquisition
The MRI data were acquired at a MRIdian MR-linac using the vendor's torso receiver coils.The treatment delivery system and MRI scanner were decoupled before imaging to operate the MR-linac in quality assurance mode to allow for orthogonal cine MRI acquisition and modification of the sequence parameters.
The following MRI data were acquired: first, the phantom motion was paused at the mid-exhale phase position and a 3D-MRI dataset was acquired with a clinical balanced steady-state free precession (bSSFP) sequence (TrueFISP; sagittal slices; slice thickness: 3 mm; in-plane resolution: 1.5 × 1.5 mm 2 ; TR/TE: 3.0/ 1.27 ms; bandwidth: 604 Hz/pixel; flip angle: 60°; acquisition matrix: 360 × 310 × 144).This simulated the breath-hold scan typically done clinically.After this scan, the breathing motion was resumed and the phantom was breathing with the same motion pattern as during 4D-CT acquisition.For each injected nodule, one timeresolved orthogonal cine MRI series (TrueFISP; slice thickness: 5 mm; in-plane resolution: 3.5 × 3.5 mm 2 ; TR/TE: 2.4/1.1 ms; bandwidth: 1000 Hz/pixel; flip angle: 60°; acquisition matrix: 100 × 100; number of averages: 1) was acquired, where the intersection line of the orthogonal slices was positioned at the approximate nodule centroid position.Each orthogonal cine MRI series consisted of 600 frames (300 in sagittal and 300 in coronal orientation) acquired over a time period of 82 s, corresponding to a frame rate of 7.3 Hz (or 3.65 Hz per sagittal/coronal pair).

Creation of continuous time-resolved estimated synthetic CTs (tresCTs)
For each of the ten gelatin nodules of the three datasets, a tresCT was created following the workflow sketched in figure 1.
For each dataset, all gelatin nodules, the lung, and the diaphragm were delineated on the mid-exhale phase (30%) image of the GT 4D-CT dataset (4D-CT 30% ).The 4D-CT 30% image was deformably registered to all remaining 4D-CT phase images using regularized B-splines with mean squared error as similarity metric using a multi-stage multi-resolution (4 stages) approach using the software Plastimatch (Shackleford et al 2010).The resulting DVFs were applied to the nodule structures defined on the 4D-CT 30% image to obtain the GT nodule positions in each of the ten breathing phases of the GT 4D-CT.
The 4D-CT 30% image was then rigidly registered to the mid-exhale phase 3D-MRI dataset using a research version of the treatment planning system RayStation (RaySearch Laboratories, Stockholm, Sweden; version 10.1.100.0) with mutual information as similarity metric.The final registration results were visually inspected in overlay plots by assessing the alignment of the inner and outer phantom walls and ten multimodality markers (MR PinPoint, Par Scientific A/S, Odense, Denmark) attached to the outer surface of the upper phantom shell.The optimal translation and rotation parameters output by the registration were applied to all GT 4D-CT phase images to get the CT and MRI data in the same frame of reference.
The 4D-CT 30% image was then deformably registered to the mid-exhale 3D-MRI using B-splines with mutual information as similarity metric and three resolution stages using Plastimatch (step 1 in figure 1) to create a synthetic 3D-CT image in the mid-exhale phase, referred to as the reference synthetic 3D-CT (sCT ref ).By choosing the mid-exhale phase as the reference image phase, we aimed at sampling the nodule motion in both directions towards the inhale phase, as well as towards the exhale phase, since the lung tumor positions observed during gated treatment at the MRIdian MR-linac are typically normally distributed around the position ) and the reference 3D-MRI dataset in the mid-exhale position are rigidly and deformably registered (DIR), and a lung density correction is applied to the resulting synthetic 3D-CT image.In step 2, the reference 3D-MRI dataset and an orthogonal cine MRI series serve as input data to the propagation method (PROP).The resulting estimated DVFs are applied to the sCT ref image, followed by a lung density correction, to create a tresCT.In step 3, this tresCT is compared to the measured GT 4D-CT, and the sCT ref image in geometric and dose calculation analyses.Measured CT data is shown in green, measured MRI data in red, and reconstructed synthetic CT data in blue boxes.
observed on the daily setup 3D-MRI scans (i.e. the breathing phase used for gating) (van Sörnsen de Koste et al 2018).To account for potential differences in the lung inflation state in the 3D-MRI dataset and the 4D-CT 30% image, the lung density correction described further below was applied to the resulting sCT ref image.The nodule contours were transferred from the 4D-CT 30% image to the sCT ref image by applying the corresponding deformation vector field (DVF) of this DIR to the individual structures.
Subsequently, the propagation method was applied to create one tresCT for each of the ten gelatin nodules of the three datasets.The propagation method is described in detail in Paganelli et al (2018b), Rabe et al (2021) and Rabe (2022).Thus, only a summary of the main steps is given here.
The mid-exhale 3D-MRI was defined as the reference 3D-MRI image (Rabe et al 2021) and was input to the propagation method together with the corresponding sCT ref image and the orthogonal cine MRI series (step 2 in figure 1).The first 15 slice pairs of the orthogonal cine MRI series were discarded to avoid frames for which the magnetization vector has not yet reached the steady state (Bieri andScheffler 2013, Rabe et al 2021).For each of the remaining 285 slice pairs, the respective sagittal and coronal slices were extracted from the 3D-MRI reference image and deformably registered to the acquired orthogonal cine MRI slice pair in independent DIRs in 2D using Plastimatch, using B-slines at six resolution stages with gradient magnitude as similarity metric.A binary image of the phantom cavity including the lung and diaphragm was used to mask the DIRs to focus the registration on the lung and to allow for sliding motion between the inner phantom walls and the lung tissue (Rabe et al 2021).The resulting DVFs were extrapolated to 3D and superimposed using the position-dependent weighting factors depending on the normal distances to the orthogonal slices introduced by Rabe et al (2021) (see equations (1)-( 4) in their publication).All deformation vectors outside of the lung and diaphragm structures defined on the sCT ref image were set to zero to prevent distortions of the static phantom shells due to extrapolation effects, since the DIRs were masked and focused on the lung and diaphragm.The resulting 3D-DVF was applied to the sCT ref image to output an estimated synthetic 3D-CT at the respective time point of the acquired orthogonal cine MRI slice pair.Finally, the lung density correction described below was applied to the resulting estimated synthetic 3D-CT.Additionally, the estimated DVFs output by the propagation method were applied to the nodule structures defined on the sCT ref to obtain the estimated nodule positions at each time point.
By repeating these steps for all 285 slice pairs (i.e.time points) for each of the ten acquired orthogonal cine MRI series (one for each of the ten gelatin nodules), a total of ten continuous tresCTs with a temporal resolution of 3.65 Hz (half the temporal resolution of the acquired orthogonal cine MRI series, since a pair of a sagittal and a coronal slice is taken as input data) and the voxel size of the GT 4D-CT (1.074 × 1.074 × 3 mm 3 ) were created.

Lung density correction
To account for lung density changes due to differences in air volume in the breathing states of (a) the 4D-CT 30% and the reference 3D-MRI dataset (step 1 in figure 1  is the determinant of the Jacobian of the DVF f of (a) the DIR of the 4D-CT 30% image to the reference 3D-MRI dataset or (b) of the estimated extrapolated DIR in the propagation method.

Breathing phase assignment of tresCTs
The geometric accuracy of each of the ten tresCTs and their usability for dose reconstruction was evaluated by comparing them to the respective GT 4D-CT images (step 3 in figure 1).The tresCTs were additionally compared to the sCT ref image to assess the added value of the tresCTs with respect to the current clinical scenario in which no time-resolved synthetic CT data is available for dose reconstruction.
To enable a breathing phase-specific comparison of the tresCT and GT 4D-CT images, the images at all time points of the tresCT had to be assigned to one of the ten breathing phase bins of the GT 4D-CT.For this purpose, a surrogate signal correlated to the breathing phase was derived by summing all pixels in a binary thresholded region of interest of the sagittal slice image of each acquired orthogonal cine MRI slice pair, which included parts of the moving superior boundary of the silicone diaphragm (Rabe et al 2021).The time points with the minimum value of summed pixels corresponded to the inhale phase (0% phase in the GT 4D-CT images).The remaining breathing phases were derived by dividing the time period between two inhale positions into ten equitemporal breathing phases.This way, each 3D-CT image of the tresCTs was assigned to one of the ten breathing phase bins.
The images in corresponding breathing phases of the tresCT and GT 4D-CT were compared in geometric and dose calculation analyses described in the following.

Geometrical analysis
The motion amplitudes of the nodule centroids in the GT 4D-CT dataset were measured from the inhale to the exhale phase.The nodule centroid positions were extracted for all ten breathing phases in the GT 4D-CT (GT nodule centroid positions) and all time steps of the tresCT (estimated nodule centroid positions).The tracking error was defined as the Euclidean distance between the estimated and GT nodule centroid positions.The tracking error was calculated for the nodule intersected by both orthogonal slices for all time points for all ten nodules of the three datasets.

Dose calculation evaluation
The dose calculation evaluation was performed with a research version of the ViewRay MRIdian treatment planning system (version 5.4.0.97).One treatment plan was created for each of the ten surrogate target lesions for the three datasets.For each treatment plan, the 4D-CT 30% image was used as the planning image, and the nodule structure intersected by the orthogonal slices was defined as the gross tumor volume (GTV).The planning target volume (PTV) was created by expanding the GTV with an isotropic margin of 5 mm, similar to common clinical practice for stereotactic MR-guided adaptive radiotherapy of lung tumors with the MRIdian MR-linac (Finazzi et al 2020b).A 6 MV flattening-filter-free step-and-shoot intensity-modulated radiotherapy (IMRT) plan with a planned dose of 5 × 11 Gy, prescribed to the PTV, was optimized while considering the effect of the 0.35 T magnetic field to mimic gated stereotactic treatment plans (Finazzi et al 2020b).A total of 14-16 beams with an angular spacing of 17.1°were set up, where angles for which the beam would enter from the non-nodule-bearing lung side were omitted.The plans were normalized to a 95% prescription dose (D p ) coverage of the PTV (V 100% 95%) (Finazzi et al 2020b).Following the recommendations of AAPM Task Group 101 for stereotactic body radiation therapy (Benedict et al 2010), dose calculation was performed on an isotropic 2 mm grid with a Monte Carlo algorithm with statistical uncertainty of 1%.
The final optimized plan was recalculated on all GT 4D-CT phases and ten estimated synthetic CTs, where one image per breathing phase was randomly sampled from the tresCT.One plan for each of the ten nodules was created and was recalculated for all ten breathing phases, yielding 100 dose distributions on both the GT and the estimated synthetic CTs.The treatment plan was additionally recalculated on the sCT ref image used as input to the propagation method.As during plan optimization, all treatment plan recalculations were performed on an isotropic 2 mm dose grid with a Monte Carlo algorithm with statistical uncertainty of 1% under consideration of the effect of the magnetic field.

Gamma pass rate analysis
For each of the ten breathing phases, the dose distributions on the GT 4D-CT images were compared to the ones on the randomly sampled tresCT images in the same breathing phase and additionally to the static dose distribution of the sCT ref image in global gamma analyses with a (2% of D p /2 mm) passing criterion for all voxels with a dose level above 10% of D p in the GT images.

Dose-volume histogram (DVH) parameter analysis
All structures defined on the 4D-CT 30% image (GTV, remaining nodules, lung, diaphragm) were rigidly copied to all other GT 4D-CT phase images, the randomly sampled tresCT images, and the sCT ref image.The DVH parameters were retrieved for these static structures defined on the GT 30% breathing phase to focus the analyses on the region with the highest dose values.The following DVH parameters were retrieved from all recalculated dose distributions: GTV D 98% , GTV D 50% , GTV D 2% , D 50% of the remaining nodules (excluding the GTV) for which the D 50% value in the planned dose distribution exceeded 11 Gy (20% of D p ), and lung V 20Gy .The median absolute differences between the GT DVH parameters in the inhale and exhale phase were computed for all analyzed DVH parameters to quantify their variation within a breathing cycle.The absolute relative deviations in percent between the DVH parameters for the GT 4D-CT dataset with respect to the corresponding DVH parameters for the tresCT dataset in the same breathing phase (Δ tresCT ), and the sCT ref image ( sCT ref D ) were calculated for all DVH parameters D listed above:

Statistical analysis
Pairwise comparisons of each of the analyzed DVH parameters were performed between Δ tresCT and sCT ref D for each breathing phase individually, and for all breathing phases combined, with a two-tailed Wilcoxon signedrank test.Furthermore, the gamma pass rates of the tresCT compared to GT 4D-CT analysis and the gamma pass rates of the sCT ref compared to GT 4D-CT analysis were also compared pairwise for each breathing phase separately, and additionally for all breathing phases combined, with a two-tailed Wilcoxon signed-rank test.All tests were performed with Python (version 3.6.5)using the implementation of the Wilcoxon signed rank test in the package scipy (scipy.stats.wilcoxon;version 1.5.4).A p-value <0.05 was considered to be statistically significant.

Geometrical analysis
The results of the geometric analysis for each of the ten nodules of the three datasets are listed in table 1.The volumes of the gelatin nodules ranged between 10-18 cc, with a mean value of 13.3 cc.With a breathing frequency of 12 cycles min −1 , and an acquisition time of 82 s, approximately 16 breathing cycles were recorded during each orthogonal cine MRI series acquisition.The inhale-to-exhale nodule centroid motion amplitudes, measured in the GT 4D-CT datasets, were in the range 3-16 mm, with a mean value of 8.2 mm.The largest mean motion amplitudes were observed for the nodules in dataset 2 (11.7 mm), followed by dataset 1 (7.3 mm), and dataset 3 (5.8mm).Averaged over all time points and all ten datasets, the (mean ± 1σ) tracking error for the nodules intersected by both orthogonal slices was (2.1 ± 0.8) mm.
Exemplary sagittal views of GT and tresCT images (dataset 1 with orthogonal slices intersecting nodule 2) and their differences in the inhale and exhale phase at two different sagittal slice positions are depicted in figure 2. In these images, a high agreement between the GT and tresCT images can be observed for both the slice intersecting the GTV, as well as a distant slice in the contralateral lung.In general, slightly larger image differences were observed in the inhale than in the exhale phase.

Dose calculation evaluation
Figure 3 shows axial views of the GT and tresCT images and recalculated dose distributions for an exemplary dataset (dataset 3 with orthogonal slices intersecting nodule 3) and the respective relative dose differences for the inhale and exhale phase.For the majority of voxels in figure 3, the relative dose differences are positive (i.e. the dose is overestimated for the tresCT with respect to the GT image) for the inhale phase, and negative (i.e. the dose is underestimated for tresCT with respect to the GT image) for the exhale phase.For 54% (57%) of the voxels in the slice depicted in figure 3 with a dose level of at least 25% of D p , the absolute dose difference was below 1%, and for 90% (95%) of the voxels the absolute dose difference was below 3% for the inhale (exhale) phase.The largest dose differences were observed in the high-dose region in the vicinity of the PTV boundary, with absolute dose differences of above 10%.

Gamma pass rate analysis
Figure 4 shows the gamma analysis results with a (2% of D p /2 mm) passing criterion as a function of the breathing phase for the tresCT and static sCT ref images.The median gamma pass rates were larger for the tresCT compared to the sCT ref images for all breathing phases except for the 20%, 30% (mid-exhale phase; breathing phase of the planning images and the reference images input to the propagation method), 70%, and 80% (midinhale phase) phases.In general, the more similar the images were to the mid-exhale phase image, the larger the median gamma pass rates were.Due to the asymmetric shape of the breathing curve, the anatomy in the exhale phase (50%) was more similar to the one in the mid-exhale phase (30%) than in the inhale phase (0%/100%), leading to the lowest median gamma pass rates for the latter for both the tresCT and the sCT ref images.

DVH parameter analysis
Figure 5 shows the cumulative DVHs for the GT, tresCT, and sCT ref images for dataset 2 where the orthogonal slices intersected nodule 2, for the GTV (nodule 2), nodules 1 and 3, and the lung for the inhale and exhale phase.For the exhale phase, the deviations of the cumulative DVHs of both the tresCT and sCT ref images with respect to the GT images were small for all analyzed structures.For the inhale phase, larger deviations for the GTV and the distant nodule 3 for sCT ref and for the GTV for tresCT with respect to the GT cumulative DVHs can be observed.No differences between the cumulative DVHs of the lung for the GT, tresCT, and sCT ref images are visible in figure 5 for neither the inhale nor the exhale phase.
The median absolute differences between the GT DVH parameters in the inhale and exhale phase, averaged over all ten datasets, were 4.1 Gy (7.4% relative difference) for the GTV D 98% , 2.9 Gy (4.5%) for the GTV D 50% , 1.7 Gy (2.5%) for the GTV D 2% , 0.3 Gy (1.5%) for the nodules D 50% , and 8.1 cc (2.0%) for the lung V 20Gy .Figure 6 shows Δ tresCT and sCT ref D for all analyzed DVH parameters and datasets as a function of the breathing phase.In general, the median Δ tresCT and sCT ref D values were largest (i.e. the largest deviations with respect to the GT values) for the breathing phases close to the inhale phase, followed by the exhale phase.The largest median Δ tresCT values were observed for the GTV D 98% with a maximum value of 2.7% for the 0% phase.For the GTV D 50% , GTV D 2% , the nodules D 50% , and the lung V 20Gy the median Δ tresCT was smaller than 2.3% for all breathing phases.
The breathing phase-specific Wilcoxon signed-rank test analysis indicated that Δ tresCT was significantly (p < 0.05) smaller than sCT ref D for the 0% and 10% phases for the GTV D 98% , and for the 0%, 10%, and 90% phases for the GTV D 50% .sCT ref D was significantly (p < 0.05) smaller than Δ tresCT only for the 30% phase for the The 30% (mid-exhale) phase is indicated as a black vertical line, representing the reference phase for the propagation method.An asterisk symbol ( * ) above the graph indicates that the gamma pass rates were significantly (p < 0.05) larger for the specific breathing phase for the tresCT with respect to the sCT ref images, and vice versa for the hash symbol (#).
lung V 20Gy .No significant differences between Δ tresCT and sCT ref D were observed for the GTV D 2% and nodules D 50% .

Discussion
The presented tresCT generation method was based on previous work by Paganelli et al (2018b) and Rabe et al (2021).The propagation method was adapted and extended to additionally input a 3D-CT to output continuous tresCTs with a temporal resolution of 3.65 Hz and the voxel size of the acquired 4D-CT (voxel size: 1.074 × 1.074 × 3 mm 3 ).The proposed method was experimentally validated with a porcine lung phantom by comparing the ten generated tresCTs to measured respiratory-correlated GT 4D-CT datasets in geometric and dose calculation analyses.
The geometrical analysis yielded a mean nodule motion amplitude of 8.2 mm.This is slightly larger than the 6 mm maximal motion that is allowed within a typically used gating window of 3 mm at the MRIdian MR-linac (van Sörnsen de Koste et al 2018, Regnery et al 2022).The average tracking error was smaller than the in-plane resolution of the cine MRI series images (3.5 × 3.5 mm 2 ) for all nodules except for the nodule with the largest motion amplitude (dataset 2, nodule 3; motion amplitude: 15.9 mm; tracking error: (3.9 ± 1.0) mm).Overall, the (mean ± 1σ) GTV tracking error averaged over all datasets and breathing phases was (2.1 ± 0.8) mm.This is slightly larger than the GTV tracking error reported by Rabe et al (2021) (median/95th percentile value of 1.5 mm/3.8 mm).The larger value found in the present study could be due to the fact that the GT and tresCT data were derived from imaging data acquired at two different scanners (CT scanner and MR-linac) instead of just the MRI scanner unit of the MRIdian, as in Rabe et al (2021).In agreement with the results described by Rabe et al (2021), the tracking error averaged over all breathing phases reported in this study tended to be larger for nodules with larger motion amplitudes (table 1).
A high agreement of the dose distributions calculated on the GT 4D-CT and randomly sampled tresCT images was found, with a median gamma pass rate of 98.6% with a (2% of D p /2 mm) passing criterion.Overall, the accuracy of the tresCT DVH parameters was high for all breathing phases with a maximum median Δ tresCT value of 2.7% for the GTV D 98% (inhale phase) and smaller than 2.3% for all other analyzed DVH parameters.The largest absolute relative deviations of up to ±10% between the dose distributions were observed in the vicinity of the PTV margin in figure 3.In this area, the dose levels are high, the dose gradients are the steepest, and small geometric offsets of the GTV can lead to pronounced differences in the local dose distributions due to the large density differences of the lung tissue and gelatin nodule.
The gamma pass rate and DVH parameter accuracy of the tresCTs with respect to the GT datasets were breathing phase-dependent.The highest accuracy was achieved for images close to the breathing phase of the reference image used in the propagation method (mid-exhale phase).In general, the dose in the lung tissue surrounding the GTV was overestimated in the inhale phase, and underestimated in the exhale phase (figure 3).A possible explanation for this is that the applied lung density correction was not sufficiently accounting for the lung density changes within a breathing cycle.
Averaged over all datasets and breathing phases, the tresCTs significantly (p < 0.05) outperformed the static sCT ref datasets in the dose calculation analyses for the gamma pass rate analysis (figure 4) and the DVH parameter GTV D 50% .The mean Δ tresCT was smaller by 25% compared to the mean sCT ref D for the GTV D 98% and approximately the same for the GTV D 2% , but these differences were not statistically significant (p > 0.05).The results suggest that even for small motion amplitudes, like the residual motion within the gating window in breath-hold-gated treatments, the tresCTs can significantly improve the dose reconstruction accuracy to the GTV with respect to today's standard method where no motion and a static image (sCT ref ) is assumed.In contrast to patients, the porcine lung phantom allows the acquisition of GT datasets for validation of experimental methods, as previously demonstrated at the MRIdian system (Rabe et al 2021), and further imaging modalities (Meijers et al 2020, Schmitz et al 2021, Bondesson et al 2022).Therefore, the measurements presented in this study are an important step towards a potential clinical implementation but are subject to a few limitations.This includes properties of the phantom, such as the lack of a heart and lung perfusion, the fact that the lung motion is purely driven by an artificial silicon diaphragm without any chest movement, and the lower lung density compared to patients.For patients, a lower density difference between the lesions and the lung tissue is expected, which potentially results in smaller dose distribution perturbations due to the magnetic field in and close to the target compared to our results observed for the porcine lung.Furthermore, the geometric errors caused by potential slightly different motion patterns at the CT scanner and MR-linac could directly impact and bias the dose calculation results.By stabilizing the lung position with an uninterruptible power supply during transport, reproducible lung tumor motion at the CT scanner and MR-linac could be achieved.The high agreement between the tresCT and GT images and dose distributions suggest that these geometrical errors were small compared to the average nodule motion amplitudes.
Furthermore, the tresCT generation method is subject to the same limitations of the propagation method, as already identified and discussed by Paganelli et al (2018b) and Rabe et al (2021).Briefly summarized, this includes DIR uncertainties, the 3D DVF extrapolation, geometric distortions of the acquired orthogonal slices, the limited spatial resolution of the currently available cine MRI sequence, and the potential localization errors related to large out-of-plane motions.Additionally, the proposed tresCT generation method requires a synthetic 3D-CT image created with a DIR of a pre-treatment 3D-CT to a 3D-MRI.This workflow step introduces additional DIR uncertainties.To minimize these, the 3D-CT and 3D-MRI need to be acquired in a similar breathing phase.This was done accordingly in the presented experiments and is standard clinical practice at the MRIdian MR-linac today as part of the routine synthetic 3D-CT generation during treatment planning.
The tresCT generation method is intended to be applied during beam-on imaging at the MRIdian MR-linac.It only requires pre-treatment 3D-CT and 3D-MRI datasets and time-resolved orthogonal cine MRI series acquired during beam delivery as input data.The first clinically implemented version of the MRIdian MR-linac only allowed cine MRI acquisition for target motion monitoring and gated beam delivery of a single sagittal slice at 4 Hz (Klüter 2019) and in a later version at 8 Hz.The latest version (called the MRIdian A3i) now allows multiplanar target tracking with cine MRI acquistion of up to three orthogonal planes (Snyder et al 2023).Thus, all of the input data required for the tresCT generation method can already be acquired in clinical routine at the MRIdian A3i MR-linac today.Therefore, the method can potentially be clinically implemented with only a few minor changes to today's workflow.
While this study was focused on the ViewRay MRIdian MR-linac, the tresCT method could potentially also be employed at the 1.5 T Elekta Unity MR-linac (Winkel et al 2019).At the Unity MR-linac, orthogonal cine MRI can also be acquired during patient treatment (Jassar et al 2023), and intrafractional dose reconstruction based on imaging data acquired during patient treatment and linac log-file data has been demonstrated in the literature (Menten et al 2020).
For both the ViewRay MRIdian and the Elekta Unity MR-linac, one research goal is to develop real-time 3D cine MRI (also named real-time 4D-MRI) sequences that can be acquired during beam delivery (Paganelli et al 2018a, Stemkens et al 2018, Kurz et al 2020, Rabe et al 2020).Such a real-time 4D-MRI dataset could be used instead of the orthogonal cine MRI data used in the tresCT generation method to derive the DVFs to warp a reference synthetic 3D-CT to create time-resolved synthetic 4D-CT datasets (Johnstone et al 2018).Grimbergen et al recently reported on a reconstruction of the delivered dose at the Unity MR-linac based on a real-time 4D-MRI (Grimbergen et al 2023).However, these sequences cannot be easily adapted and employed at the MRIdian MR-linac with its low magnetic field strength of 0.35 T and restricted MRI hardware (Stemkens et al 2018, Klüter 2019).Therefore, the presented tresCT generation method might represent a valuable intermediate method for intrafractional reconstruction of the delivered dose for moving lung tumors until such sequences become potentially clinically available at the MRIdian MR-linac in the future.

Conclusions
We proposed a method to create continuous time-resolved estimated synthetic 4D-CTs (tresCTs) using 3D-CT, 3D-MRI, and orthogonal cine MRI series as input data.The method was experimentally validated with a porcine lung phantom.The resulting tresCTs achieved high agreement with measured GT respiratory-correlated 4D-CT datasets in terms of geometric and dose calculation accuracy.The proposed tresCT generation method could be a valuable tool for retrospective dose reconstruction or online intrafractional dose accumulation during MRgRT of lung cancer patients treated at the ViewRay MRIdian MR-linac.The simplicity of the method and the similarity of the proposed workflow to today's clinical workflow could potentially enable a low-threshold clinical implementation in the future.

Figure 1 .
Figure1.Workflow of tresCT creation and comparison.In step 1, the mid-exhale phase of the measured GT 4D-CT dataset (4D-CT 30% ) and the reference 3D-MRI dataset in the mid-exhale position are rigidly and deformably registered (DIR), and a lung density correction is applied to the resulting synthetic 3D-CT image.In step 2, the reference 3D-MRI dataset and an orthogonal cine MRI series serve as input data to the propagation method (PROP).The resulting estimated DVFs are applied to the sCT ref image, followed by a lung density correction, to create a tresCT.In step 3, this tresCT is compared to the measured GT 4D-CT, and the sCT ref image in geometric and dose calculation analyses.Measured CT data is shown in green, measured MRI data in red, and reconstructed synthetic CT data in blue boxes.
) or (b) the sCT ref image and the acquired orthogonal cine MRI slice pairs (step 2 in figure 1), a scaling of the CT numbers of (a) the sCT ref image or (b) the tresCT was performed, following the method by Sarrut et al (2006) (equation (5) in their publication): the CT numbers in (a) the sCT ref image or (b) the tresCT to be corrected, HU 2 are the CT numbers in (a) the 4D-CT 30% image or (b) the sCT ref image, and ( ) det f

Figure 2 .
Figure 2. Visual comparison of GT and tresCT images.The figure shows exemplary sagittal views of the GT (upper row) and tresCT (center row) images, and their differences (bottom row).The images are shown in the inhale (first and third column) and exhale phase (second and fourth column) at two different sagittal slice positions (a) and (b).The position of the acquired coronal slice used as input to the propagation method is indicated as a yellow vertical line.All contours are shown in the mid-exhale phase, as defined on the 4D-CT 30% image.(a) A sagittal slice at the position of the acquired sagittal slice in the right lung, intersecting the GTV (orange contour), the PTV (red), and another nodule (green) is shown.The diaphragm is shown in blue.(b) A sagittal slice located in the left lung at a distance of 14.1 cm to the acquired sagittal slice and intersecting a distant nodule (cyan) is shown.

Figure 3 .
Figure 3. Recalculated dose distributions on GT and tresCT.The figure shows exemplary axial views of the images and dose distributions of the GT dataset (left column), the tresCT (center column), and the corresponding relative dose differences (right column) for the inhale (top row) and exhale phase (bottom row).The dose is plotted relative to the prescribed dose D p = 55 Gy.The dose differences are plotted relative to the prescribed dose: (D tresCT −D GT )/D p .Doses below 25% of D p and dose differences below an absolute value of 1% were masked for clearer visualization.The PTV is shown as a red contour and the diaphragm as a blue contour (contours as defined on the 4D-CT 30% image).

Figure 4 .
Figure4.Gamma pass rate analysis results.The solid lines indicate the median gamma pass rates for the dose distributions on the tresCT (blue) and sCT ref (red) images with respect to the GT dose distributions as a function of the breathing phase for all ten treatment plans (ten different tracked nodules).The shaded areas mark the minimum and maximum values for each breathing phase.The 30% (mid-exhale) phase is indicated as a black vertical line, representing the reference phase for the propagation method.An asterisk symbol ( * ) above the graph indicates that the gamma pass rates were significantly (p < 0.05) larger for the specific breathing phase for the tresCT with respect to the sCT ref images, and vice versa for the hash symbol (#).

Figure 5 .
Figure 5. Cumulative DVHs of analyzed structures.The cumulative DVHs for the analyzed structures are plotted for the GT (solid line), tresCT (dashed line), and sCT ref (dotted line) images for the inhale (left) and exhale (right) phases.

Figure 6 .
Figure 6.DVH parameter analysis results.The median Δ tresCT and sCT ref Daveraged over all ten datasets are plotted as solid lines as a function of the breathing phase for all analyzed DVH parameter (subplots a-e).The shaded areas indicate the minimum and maximum values for each breathing phase.The 30% (mid-exhale) phase is indicated as a black vertical line, representing the reference phase for the propagation method.An asterisk symbol ( * ) above the graph indicates that Δ tresCT was significantly (p < 0.05) smaller than sCT ref D for the specific breathing phase, and vice versa for the hash symbol (#).

Table 1 .
Geometrical analysis results.The GTV size, motion amplitude A motion , and tracking error are listed for each of the ten nodules of the three datasets individually and averaged over all nodules.
The mean Δ tresCT and sCT ref D values for all DVH parameters, averaged over all breathing phases and datasets are summarized in table 2. The mean Δ tresCT values were smaller by 25% and 40% compared to the sCT ref D values for the GTV D 98% and GTV D 50% .The difference was statistically significant (p < 0.01) for the GTV D 50% .For the GTV D 2% and nodules D 50% , the mean Δ tresCT was slightly larger than the sCT ref D value, but the difference was not statistically

Table 2 .
Mean absolute relative deviations of DVH parameters.The mean Δ tresCT and sCT refDvalues and their differences, averaged over all breathing phases and datasets are reported in %.The p-values of the Wilcoxon signed rank test are reported in the second to last column.An asterisk symbol ( * ) indicates that the deviations were significantly (p < 0.05) smaller for the the tresCT with respect to the sCT ref images, and vice versa for the hash symbol (#).