Marc Derieppe
email: marc.derieppe@gmail.com

Clemens Bos

M De Greef

Chrit Moonen

Baudouin Denis De Senneville

Correction for Photobleaching in Dynamic Fluorescence Microscopy: Application in the Assessment of Pharmacokinetic Parameters in Ultrasound-Mediated Drug Delivery

Keywords: drug delivery, pharmacokinetic parameters, photobleaching, ultrasound bioeffects, dynamic fluorescence microscopy

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Local drug delivery in oncology aims at increasing the dose of an anticancer agent at the tumour site while limiting its concentration in the general circulation [START_REF] Allen | Drug delivery systems: entering the mainstream[END_REF][START_REF] Yeh | Cardiovascular Complications of Cancer Therapy: Diagnosis[END_REF]. Efficient delivery remains challenging owing to the existence of biological barriers, such as the endothelial barrier and tumor-cell plasma membranes for agents that need to reach the intracellular compartment to induce its cytotoxic effect.

To increase the permeability of these biological barriers towards compounds exhibiting a low lipophilicity, strategies exploiting acoustic cavitation generated by ultrasound waves (US) have been devised in a localized and non-ionizing fashion [3][START_REF] Miller | A review of in vitro bioeffects of inertial ultrasonic cavitation from a mechanistic perspective[END_REF][START_REF] Mitragotri | Healing sound: the use of ultrasound in drug delivery and other therapeutic applications[END_REF][START_REF] Haar | Therapeutic applications of ultrasound[END_REF][START_REF] Haar | Ultrasound mediated drug delivery: A 21st century phoenix?[END_REF]. Injected intravenously, these gas-filled micrometric particles remain located in the vasculature and induce permeability of the endothelial cells, thus facilitating drug extravasation of therapeutic agents that exhibit poor transendothelial transport. Monitoring of these processes in vitro at the cellular level and in vivo at the endothelialbarrier level is of special interest to optimize US-mediated delivery protocols and gain insight in US bioeffects.

For this purpose, fluorescence microscopy provides a spatial resolution suitable for cell imaging, and a wide range of fluorescent probes is available to observe mechanisms selectively. As an example, the use of cell-impermeable fluorophores demonstrating a fluorescence signal enhancement when bound to DNA was reported; these fluorophores serve as "smart sensors" for the investigation of US-and microbubble-mediated plasma membrane permeabilization [START_REF] Deckers | A fluorescent chromophore TOTO-3 as a 'smart probe' for the assessment of ultrasound-mediated local drug delivery in vivo[END_REF]. With a molecular weight close to that of anticancer agents, these fluorophores, e.g. SYTOX Green, 600 Da or TO-PRO-3, 671 Da, can be considered models for low molecular weight anticancer drugs.

Capable of collecting images with a high temporal resolution, modern fluorescence microscopy systems paved the way for a broad spectrum of applications providing insights in the spatiotemporal dynamics of biological systems [START_REF] Fricker | Quantitative Fluorescence Microscopy: From Art to[END_REF][START_REF] Muzzey | Quantitative Time-Lapse Fluorescence Microscopy in Single Cells[END_REF]. For instance, fluorescence imaging of molecular motors in the cytoplasm [START_REF] Stepanova | Visualization of Microtubule Growth in Cultured Neurons via the Use of EB3-GFP (End-Binding Protein 3-Green Fluorescent Protein)[END_REF][START_REF] Jaqaman | Robust single-particle tracking in live-cell time-lapse sequences[END_REF] allowed collecting speed and trajectory estimates, and could relate them to biological conditions. Unfortunately, direct investigation of the fluorescence signal in intensity-based fluorescence microscopy remains challenging due to photobleaching (PB).

PB stems from the photon-induced chemical destruction of a fluorescent molecule that is subjected to excitation light, resulting in the irreversible reduction of fluorescence signal intensity [START_REF] Rosenthal | Optical Imaging of Cancer: Clinical Applications[END_REF]. Well-described in homogeneous solutions with a controlled chemical environment [START_REF] Song | Photobleaching kinetics of fluorescein in quantitative fluorescence microscopy[END_REF], the PB process remains difficult to assess in biological systems where fluorophores are exposed to a multimolecular chemical microenvironment. This limits systematic quantitative analysis of fluorescence signal intensity data. Even if a number of studies showed that PB is not a fully monoexponential process [START_REF] Song | Photobleaching kinetics of fluorescein in quantitative fluorescence microscopy[END_REF][START_REF] Benson | Digital imaging fluorescence microscopy: spatial heterogeneity of photobleaching rate constants in individual cells[END_REF][START_REF] Rigaut | High-resolution three-dimensional images from confocal scanning laser microscopy. Quantitative study and mathematical correction of the effects from bleaching and fluorescence attenuation in depth Anal Quant[END_REF], it is commonly assumed that a single-exponential process describes PB in good approximation [START_REF] Smal | Advanced particle filtering for multiple object tracking in dynamic fluorescence microscopy images 4th[END_REF].

In our previous studies, we demonstrated in vitro the feasibility to monitor in real-time the USand microbubble-mediated permeabilization of a SYTOX Green model drug using a fibered confocal fluorescence microscopy (FCFM) system that avoids the geometrical constraints of standard microscope setups. Having real-time data available should permit to characterize the model drug uptake quantitatively in terms of the uptake kinetics [START_REF] Derieppe | Real-Time Assessment of Ultrasound-Mediated Drug Delivery Using Fibered Confocal Fluorescence[END_REF]. Thus, we implemented an automated postprocessing pipeline to assess the pharmacokinetic parameters of this model drug, including the detection of cells, cell-tracking over the full sequence, and statistical analyses of the cell population in the field of view [START_REF] Derieppe | Tracking of Cell Nuclei for Assessment of In Vitro Uptake Kinetics in Ultrasound-Mediated Drug Delivery Using Fibered Confocal Fluorescence[END_REF]. However, assessment of pharmacokinetic parameters was hindered by PB inherent of dynamic fluorescence microscopy. In this article, we assess the impact of PB on the observed pharmacokinetic parameters, and evaluate the use of three-compartment models, that include a rate constant to reflect PB.

Materials and Methods

Setup and Experimental Protocol

In this study, dynamic fluorescence data during and after sonoporation were acquired using the protocol presented in Derieppe et al. [START_REF] Derieppe | Real-Time Assessment of Ultrasound-Mediated Drug Delivery Using Fibered Confocal Fluorescence[END_REF]. Briefly, the in vitro setup comprised a 37°C water bath containing a 1.5 MHz unfocused mono-element transducer (diameter 20 mm) with the following US settings: 0.9 MPa peak-to-peak pressure, 10 % duty cycle, 1 kHz pulse repetition frequency, 5 s exposure time. These US settings do not induce any temperature rise.

Three days prior to the experiment, 500,000 C6 rat glioma cells were seeded in a 10 mL cell culture chamber (OptiCell TM , Thermo Fischer Scientific, Rochester, NY, USA) that is ultrasound and optically compatible. On the day of experiment, cells were at 90 % confluence. A commercially available ultrasound contrast agent (SonoVue ® , Bracco, Milan, Italy) was added to the cells at a 20 microbubble-to-cell ratio; this microbubble density was chosen to optimize the uptake yield, i.e., the highest percentage of cells taking up the dye, while keeping a high cell viability (data not shown). Ten seconds after starting the recording, the cell monolayer was exposed to US waves for 5 s. Making contact with the upper wall of the culture chamber, the tip of the FCFM microprobe (Z-probe -Cellvizio ® Dual Band, Mauna Kea Technologies, Paris, France) allowed monitoring the uptake of SYTOX Green fluorescent dye (excitation/emission = 504/523 nm, molecular weight = 600 Da) by means of the fluorescence signal intensity. To evaluate the impact of PB, the experiment was conducted using 0.5 mW, 1.0 mW and 1.6 mW of laser power, which are 25 %, 50 % and 75 % of the maximum laser power, respectively. The lower light intensity was chosen to detect cells first on lower signal-to-noise-ratio images at the early stage of the ultrasound-and microbubble-mediated cell uptake.

Signal Acquisition and Cell Tracking

To monitor US-and microbubble-mediated cell uptake, the fluorescence signal time curve for each individual cell was obtained at a frame rate of 8.5 Hz, i.e. the lowest allowed by the FCFM hardware, as described in Derieppe et al. [START_REF] Derieppe | Tracking of Cell Nuclei for Assessment of In Vitro Uptake Kinetics in Ultrasound-Mediated Drug Delivery Using Fibered Confocal Fluorescence[END_REF]. Image analysis was conducted using MATLAB ® 2013 (MathWorks, USA). Briefly, a first series of 6 minutes (min) was acquired during which cells were sonicated for 5 s, ten seconds after the acquisition was started (figure 1a). In these images, detection of the cell nuclei was performed using a non-local means approach followed by an intensity threshold for nucleus delineation. The nuclei were then labelled in each frame and automatically tracked over the time series using the iterative closest point method [START_REF] Cressie | Statistics for Spatial Data Wiley-Interscience Revised Edition[END_REF] applied between successive images, which is justified since the acquisition frame rate is high compared to the cell motion. Although the fluorescence signal is only collected in the cell nucleus, in the following, it will be thus referred to as the "cell signal". Evaluation of PB was carried out for the same field of view by acquiring a second sequence after the fluorescence signal further equilibrated (figure 1b). Pilot experiments allowed measuring uptake time constants in the range of 3 minutes (two-compartment model). Assuming that the equilibration is reached five times later than the characteristic uptake time constants, equilibration for most cells was sufficiently achieved 15 minutes after US exposure (figure 1c) when a second series was started in which the fluorescence signal decay observed originates from PB (figure 1b).

Mathematical Models and Computation of the Kinetic Parameters

As suggested before [START_REF] Derieppe | Real-Time Assessment of Ultrasound-Mediated Drug Delivery Using Fibered Confocal Fluorescence[END_REF], the uptake of the model drug in the first image series can be described as a two-compartment model with an extracellular compartment, which initially contains all SYTOX Green molecules, and an intracellular compartment, taking up the model drug during and after US-and microbubble-mediated sonopermeabilization. After binding to DNA, SYTOX Green fluorescence intensity increases by two orders of magnitude, thus reflecting cell uptake of SYTOX Green. Thus, for the two-compartment model, the time dependence of the fluorescence signal intensity is as follows:

(1) where is the asymptotic signal enhancement, T the time of fluorescence signal onset, and k up the uptake rate.

In the second image series, PB was evaluated by fitting the fluorescence intensity profiles with a mono-exponential decay, as suggested in the literature [START_REF] Smal | Advanced particle filtering for multiple object tracking in dynamic fluorescence microscopy images 4th[END_REF][START_REF] Song | Influence of fluorochrome labeling density on the photobleaching kinetics of fluorescein in microscopy[END_REF]:

(

) 2
where I 0 is the fluorescence signal at the onset of the acquisition, and k pb the PB rate. This model will be named "PB model" in this study.

To take PB into account during the uptake phase when fluorescence signal increases, a threecompartment model can be used where we add an efflux rate of fluorophore molecules subjected to the irreversible damage of PB. The time dependence of the fluorescence signal intensity is then given by:

(

where is the asymptotic fluorescence signal, T the time of fluorescence signal onset, k up the uptake rate, and k pb the PB rate.

The uptake rates and the PB rates were computed using a Levenberg-Marquardt least-squares fit applied to the fluorescence intensity in each individual cell. The model was considered a good description of the enhancement when the goodness of the fit as evaluated by the Pearson's correlation coefficient (r 2 ) was greater than 0.95.

Abbreviated 2CM in this study, the application of the two-compartment model consisted of fitting the uptake signal profiles of the individual cells with , T and k up as free fit parameters using the signal equation as given by Equation 1. The three-compartment model (Eq. 3) was applied on the uptake signal profile by considering two approaches. First, using the PB model on the data from the equilibrium phase, the PB rate was evaluated in each individual cell experimentally (figure 1b). For each cell, this PB rate then served as an input for a fit of the data to the three-compartment model (Eq.

3), keeping it fixed, while including , T and k up as fit parameters. Abbreviated 3CM ref , this model is the reference method in this study. Second, the three-compartment model was used to fit the data of the uptake phase directly, with all the parameters free in the fit (3CM all ).

Influence of the Acquisition Time -Sensitivity of the Three-compartment Model

Sensitivity of the three-compartment model was evaluated, in particular the influence of errors in determining k pb on the estimation of k up . The estimated uptake rate was noted . To this end, a synthetic uptake signal profile was generated with an asymptotic signal set to 1000 arbitrary units (a.u.), a signal onset at 10 s, a chosen uptake rate k up and a PB rate k pb . Then, the sensitivity of , as determined by fitting of the signal profile with 3CM, to deviations from the chosen PB rate was investigated. This was evaluated by each time generating a randomly perturbed PB rate from a

Gaussian distribution centered at the true value of k pb , with a standard deviation of 25 % k pb . The sensitivity was evaluated in a range of PB rates (at 25 %, 50 % and 75 % of laser power) and uptake rates found experimentally, i.e. from 0.1 10 -3 s -1 to 3.5 10 -3 s -1 , and from 4.0 10 -3 s -1 to 50.0 10 -3 s -1 , respectively. The simulation was performed using 3 min, 6 min and 9 min of acquisition time for the collection of the uptake profile. For each set of kinetic parameters, the model sensitivity was estimated by assessing the normalized error, defined as the difference of the median value of the estimates and the chosen k up , divided by the chosen k up , and the normalized precision, i.e. the interquartile range of divided by k up .

Comparison of the methods

For each of the three laser powers, the uptake rate estimates obtained from the 2CM (Eq. 1) and the 3CM all (Eq. 1) were compared to the 3CM ref by using a Bland-Altman plot representing the difference of the same estimate measured with 2 different methods as a function of their average value [START_REF] Bland | Statistical methods for assessing agreement between two methods of clinical measurement[END_REF][23].

Statistical Analysis

Data are presented as median (interquartile range). Statistics have been conducted using GraphPad

Prism ® software (San Diego, CA, USA). PB rates evaluated experimentally at 25 %, 50 % and 75 % of laser power were compared one-by-one using the unpaired non-parametric Mann-Whitney test (MW) (figure 3). Uptake rates as well as the PB rates obtained using the different models at different laser powers were compared using Wilcoxon's paired signed rank test (WX). The results were considered significant when the p-value (p) was lower than 0.05.

Results

Real-Time Imaging of Model Drug Uptake and Experimental Evaluation of Photobleaching

The first 6-minute acquisition shows that exposure of the cells to US waves leads to a fluorescence signal enhancement that starts at the onset of US, and continues after the end of US exposure (figure 1a, Fig. 2a-d), as was described previously by Derieppe et. al. 2013 [START_REF] Derieppe | Real-Time Assessment of Ultrasound-Mediated Drug Delivery Using Fibered Confocal Fluorescence[END_REF]. Then, after the SYTOX Green model drug had further equilibrated, the PB profile was acquired in the same field of view (figure 1b, Fig. 2e-h). The fluorescence signal at 15 minutes (figure 1b) was higher than that at 6 minutes (figure 1a), owing to influx of fluorophores in the cell nuclei binding to DNA, spared from PB between 6 minutes and 15 minutes. Using the PB model (Eq. 2), the PB rates collected cell-by-cell significantly increased (MW, p < 0.05) as a function of the laser power: 1.0 10 -3 s -1 (0.4 10 -3 s -1 ), 2.1 10 -3 s -1 (0.7 10 -3 s -1 ) and 2.9 10 -3 s - 1 (0.9 10 -3 s -1 ) at 25 %, 50 % and 75 % of laser power, respectively, indicating a clear impact of the laser power on SYTOX Green PB (figure 3).

Sensitivity of the Three-compartment Model and Influence of the Acquisition Time

In the range of uptake rates and PB rates obtained experimentally using the 3CM ref (Eq. 3), the sensitivity of the three-compartment model was evaluated for 3 time windows. The normalized errors of were small, and showed no clear dependence on the true values of k up and k pb . The maximum normalized error found was 0.6 %, 0.2 % and 0.2 % for acquisitions of 3 minutes, 6 minutes and 9 minutes, respectively, indicating that for uptake rates in our experiments it is beneficial to acquire data for at least 6 min, but that the normalized error is not lowered when the acquisition time is further to 9 min (figure 4a,c,e). The precision of was better than 8.2 %, 3.3 % and 2.7 %, for all simulated k up and k pb (figure 4b,d,f), for acquisition windows of 3 minutes, 6 minutes and 9 minutes, respectively. For the short acquisition window, the variability in the estimated uptake rates was notably higher for low  ˆ k up uptake rates, and the more so when PB rate was low. In addition, for all observation durations, variability increased for combinations of high PB rate and low uptake rate, i.e. when k up  k pb . If, for our purposes in our experiment, a precision of 5 % is regarded as satisfactory, it is thus sufficient to acquire data using a 6-minute time window.

Comparison of the methods

After the uptake and PB profiles for individual cells were extracted from the data, results from the 2CM and the 3CM all methods were compared to the 3CM ref method, which is our reference method in this study (see section "Real-Time Imaging of Model Drug Uptake and Experimental Evaluation of Photobleaching") (figure 5a-c). The Bland-Altman plot indicates a proportional bias with an increasing overestimation of the uptake rate by the 2CM when the average value increases. The overestimation increased with laser power and the PB induced, as shown by the stronger slope of the line of best fit: 0.33, 0.66 and 1.01 at 25 %, 50 % and 75 % of laser power, respectively (figure 5a-c, dark line). In addition, the scattergram expands when the average value increases, which is characteristic of heteroscedasticity. This spread increases linearly when the uptake rate increases (figure 5a-c, red lines).

For the majority of the points in the low range 0 -10 10 -3 s -1 , comparison of the 3CM all with the 3CM ref indicates a proportional bias with a negative slope, indicating that the 3CM all underestimates the uptake rate (figure S1a-c). For the other data points with an average uptake rate larger than 10.0 10 -3 s -1 , no clear trend can be identified. Considering all the points, difference values display an average of -1.7 10 -3 s -1 , -0.5 10 -3 s -1 , and 2.8 10 -3 s -1 at 25 %, 50 % and 75 % of laser power, respectively (figure S1a-c, black horizontal line).

Distribution of Pharmacokinetic Parameters

First, the uptake rates calculated by the 3CM ref at 25 %, 50 % and 75 % of laser power were in agreement (MW, p > 0.05), indicating that the 3-compartment model allowed assessing the uptake rates of the cells irrespective of the amount of PB (figure 6a-c and Table S1). Conversely, the 2CM that was commonly used before this study showed significant differences (MW, *p < 0.05) in the uptake rates obtained at different laser powers. The median values doubled and the interquartile ranges increased by a factor of 2.5 from 25 % to 75 % of laser power, i.e.

with a stronger PB (figure 6a-c), suggesting a limited accuracy of the 2CM.

With the 3CM all , which is based on the 3-compartment model (Eq. 3) and does not have experimentally obtained independent knowledge on the PB rate, the uptake rate estimates at 25 % of laser power, were significantly lower than those at 50 % (MW, *p < 0.05), however, no difference was observed in the uptake rate estimates between 50 % and 75 % of laser power (figure 6a-c and Table S1). As also observed for the 2CM, the interquartile range increases when the PB is higher. Moreover, a significant difference is observed between the 3CM ref and the 3CM all at 25 % of laser power (WX, *p < 0.05) (figure 6a), but similar uptake rate estimates were found with the 3CM all compared to the 3CM ref (WX, p > 0.05) at 50 % and 75 % of laser power (figure 6b-c and Table S1). This suggests that the 3CM all becomes more robust when exciting the fluorophores with a significant laser power. If we then compare the PB rates obtained by the 3CM all with those obtained in a separate PB measurement as used in 3CM ref , similar estimates (WX, p > 0.05), were found only at 75 % of laser power, again indicating that the PB must be high enough to be sensed in the fitting procedure (figure S2a-c).

Discussion

PB is a frequently occurring problem in applications of intensity-based dynamic fluorescence microscopy [25]. For instance, in single-particle tracking applications, signal-intensity based detection of single particles is hindered by PB [26], thus requiring its evaluation and methods to compensate for it. PB compensation is all the more important when the quantity of interest is directly derived from fluorescence signal intensity. In our specific application, PB had an impact on pharmacokinetic parameter estimates collected in US-and microbubble-mediated plasma membrane permeabilization using FCFM in vitro. A previous study proposed a two-compartment model to describe the cell uptake of SYTOX Green model drug triggered by US and microbubbles [START_REF] Derieppe | Real-Time Assessment of Ultrasound-Mediated Drug Delivery Using Fibered Confocal Fluorescence[END_REF]; to model the irreversible destruction of fluorophores by PB, a third compartment taking up the fluorophores bound to DNA was added, and it was demonstrated that incorporating PB led to more consistent results, independent of the applied laser power. In the context of US-and microbubble-mediated drug delivery, this study is the first to analyze systematically PB impact and compensate for it on the uptake kinetic parameter estimates of a model drug measured by dynamic confocal fluorescence microscopy.

After the acquisition of the uptake profile and after the SYTOX Green fluorescence signal had further equilibrated, the collection of the PB profile in each individual cell was performed in the same field of view, ensuring the same optical path. The PB rates were measured larger when the laser power increased, but overall lower than 5.0 10 -3 s -1 . In addition, the heterogeneity of the PB rate distribution increased with a greater laser power. In the scope of this study, for the sake of stability of the fit and according to the PB profiles, a mono-exponential model was chosen.

Using simulations, the numerical stability of the uptake rate estimates to wrong estimations of the PB rate was found to be low. We showed that it is advantageous to set the observation window of the imaging protocol at least twice the characteristic time constant, especially when both k up and k pb are fitted from this data simultaneously.

The 2CM and the 3CM all were compared to the 3CM ref , which integrates independent knowledge by means of an experimental PB rate. The latter was the reference method of this study.

We observed that the 2CM introduces a bias, that is, linearly proportional to the uptake rate, i.e. the larger the uptake rates in a cell population, the higher their overestimation by the 2CM. This overestimation is also observed in the synthetic data exemplified in Figure S3, and can be explained by the inability for the 2CM to fit adequately the local maximum in the fluorescence signal profile, that is introduced by PB. Therefore, the least-squares fit based on the 2CM solves the minimization problem by overestimating the uptake rate. It should also be noted that this overestimation is even larger when the PB is higher, as reflected by the slope of the proportional bias as a function of laser power. These observations indicated that PB has an impact on the uptake rate estimates, and must be taken into account to compensate for it. Finally, the methods were compared for their potential to accurately assess uptake rate in a cell population. Using PB rates from a separate measurement phase, 3CM ref found uptake rates irrespective of the laser power applied, which demonstrates adequate PB correction for this three-compartment model. To investigate if it would be possible to obtain uptake and photobleaching rate from a single acquisition window, we compared the 3CM all , in which all parameters are free in the fit procedure, to the 3CM ref .

The Bland-Altman plot did not indicate a clear trend, but an increasing underestimation is noticeable for uptake rate estimates below 10.0 10 -3 s -1 for the 3 laser powers. The cause of this underestimation and the range where it is applicable remain unclear, but at 25 % and 50 % of laser power, this underestimation of the uptake rates using the 3CM all is associated with an overestimation of the PB rates compared to those found using the 3CM ref .

In addition, the scattergrams for each of the laser powers showed no pattern for uptake rates greater than 10.0 10 -3 s -1 .

When PB and uptake are fitted simultaneously from a time series, it is necessary to perform a pilot experiment and, from the resulting time constants, set the time of acquisition so the sequence is long enough and the PB rate must be high enough, but in a range of laser power where no heat is generated, to have substantial effect in the last part of the sequence and thus be adequately sensed by the 3CM all model. This finding helps us to design future imaging protocols related to this application.

Interestingly, this runs counter to the common instruction in fluorescence microscopy advising to set the intensity of the excitation light as low as possible to limit PB.

Conclusion

This study established a methodology to evaluate PB, investigate its impact on kinetic parameter estimates measured by dynamic fluorescence microscopy, and compensate for PB using a three-

Figure 1 .

 1 Figure 1. Imaging timeline (c) with a representative fluorescence signal profile of the uptake (a) and

The

  2CM was compared to the 3CM ref by determining the line of best fit with a first-order model as described in Ludbrook et. al. 2010 [24]. The heteroscedasticity was represented by V-shaped standarddeviation lines that were drawn by fitting linearly the absolute values of the residuals, i.e. the observed values by the predicted ones, against the average values [24]. The resulting coefficients for regression were adjusted by multiplying them by √(ᴨ/2). It yields increasing standarddeviation values when the average values increase.

Figure 2 .

 2 Figure 2. Snapshots of a representative monitoring of the US-and microbubble-mediated SYTOX

Figure 3 .

 3 Figure 3. Experimental evaluation of PB rates cell-by-cell (*p < 0.05). The photobleaching rate

Figure 4 .

 4 Figure 4. Sensitivity maps of the three-compartment model. a.c.e. Normalized error with a 3 (a), 6 (c)

Figure 5 .

 5 Figure 5. Bland-Altman plots of the uptake rate estimates using the 2CM compared to the 3CM ref at 25 % (a), 50 % (b) and 75 % (c) of laser power. The scatter plots show a proportional bias with

Figure 6 .

 6 Figure 6. Box plots of the uptake rates obtained with the 3CM ref , the 3CM all and the 2CM at 25 % (a), 50 % (b) and 75 % (c) of laser power. *p < 0.05. Non-significant differences (NS): p > 0.05.

  compartment model. It was demonstrated that a 3CM model helped to improve the accuracy of the uptake pharmacokinetic parameter estimates in a cell population subjected to microbubble-assisted sonopermeabilization by US waves.

Table S1 .

 S1 Table of the uptake rate estimates obtained with the methods 3CM ref , 2CM and 3CM all . median (interquartile range).

				Method
	Uptake rates (10 -3 s -1 )
			3CM ref	2CM	3CM all
		25 %	5.5 (5.1)	7.8 (7.6)	3.2 (2.1)
	Laser Power	50 %	6.5 (5.1)	12.5 (11.7)	4.3 (5.0)
		75 %	6.1 (5.9)	15.6 (19.0)	5.2 (12.1)

Table S2 .

 S2 Table of the photobleaching rate estimates obtained with the methods 3CM ref and 3CM all . median (interquartile range).

	Photobleaching rates (

10 -3 s -1 ) Method 3CM ref 3CM all Laser Power 25 % 1.0 (0.4) 2.0 (1.9) 50 % 2.1 (0.7) 2.3 (2.0)

	75 %	2.9 (0.9)	2.8 (2.2)

Acknowledgments

This study was funded by the ERC project 268906 "Sound Pharma" (Prof. C. Moonen). The authors thank Dr. Sébastien Benzekry (IMB, UMR 5251, CNRS, Université de Bordeaux, Bordeaux, France) for the fruitful discussion on the theoretical robustness of the 3-compartment model.

Figure S3. Examples of uptake rates obtained with the 2CM from synthetic signal profiles generated with 5.0 10 -3 s -1 (a), 20.0 10 -3 s -1 (b) and 40.0 10 -3 s -1 (c) of uptake rates, a photobleaching rate of 1.0 10 - 3 s -1 (corresponding to 25% of laser power), and an asymptote of 2000 a.u. applied in each case.


Markus Glitzner
email: m.glitzner@umcutrecht.nl

Baudouin Denis De Senneville

Jan J W Lagendijk

Bas W Raaymakers

Sjoerd P M Crijns

On-line 3D motion estimation using low resolution MRI

Keywords: deformable image registration, local deformations, resampling, data sufficiency Phys. Med. Biol

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Recent developments enabled on-line magnetic resonance imaging (MRI) during radiotherapy [START_REF] Raaymakers | Integrating a 1.5 T MRI scanner with a 6 MV accelerator: proof of concept[END_REF][START_REF] Lagendijk | The Magnetic Resonance Imaging-Linac System[END_REF]. In further consequence, new possibilities for elements in the feedback chain of adaptive radiotherapy (ART) arise: in contrast to currently available real-time positioning devices which provide only surrogate displacement data, MRI returns 3D anatomical information of both target volume and organs at risk (OAR).

After acquisition, 3D displacements need to be estimated from the dynamic volumes using image-based motion estimation algorithms [START_REF] Oliveira | Medical image registration: a review[END_REF]. The displacements become 3D deformation vector field (DVF) which can be used as feedback variable for e.g. dose calculations or target tracking during intervention (de Senneville et al. 2011). However, MRI is an inherently slow image modality and acquisition speed depends directly on image resolution. Additionally, MRI acquisition speedup is limited by both technical (gradient strength/noise bandwidth) and physiological (specific absorption rate (SAR), peripheral nerve stimulation (PNS)) constraints (?). Similarly, model based acceleration techniques such as parallel imaging [START_REF] Blaimer | GRAPPA: how to choose the optimal method[END_REF]) or compressed sensing [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF]) are limited by the validity of their underlying model and the demand for reconstruction time.

The aspiration of this work is to investigate the impact of image resolution on the quality of motion estimation. Potentially, with lower sufficient resolution for comparable motion estimation, magnetic resonance (MR) acquisition could be accelerated significantly. The gain in acquisition speed would enable applications which demand for high-frequency, real-time MRI, such as 3D target and dose tracking for ART. To avoid image modulation from breathing motion, respiratory gating was employed with a gating window of 5mm using navigator echoes on the interface between lung and left liver-dome. Gating in end-of-expiration was chosen because of an expected low residual breathing motion uncertainty for abdominal sites [START_REF] Berbeco | Residual motion of lung tumors in end-of-inhale respiratory gated radiotherapy based on external surrogates[END_REF].

Methods

Since isotropic resolution is a prerequisite of the utilized image processing, the acquisition was tailored to acquire data with (2.5mm) 3 voxel size. k-space was filled in 3D manner in order to avoid malicious slice effects, such as slice overlap and inter-slice motion. Additionally, 3D k-space sampling yields higher signal to noise ratio (SNR) compared to sliced acquisitions.

The navigator triggered an echo planar imaging (EPI) readout with an EPI-factor of 21, T R /T E /α = 47ms/9.8ms/25 • . Binomial pulses (1 -2 -1) [START_REF] Hore | Solvent suppression in fourier transform nuclear magnetic resonance[END_REF] were used for selective water excitation, concurrently suppressing the fat signal. Furthermore, moderate sensitivity encoding (SENSE) (P = 1.5, S=1.5) was used to further accelerate imaging.

In

Image Processing

To simulate coarser sampling, downsampled data was synthetically generated using MATLAB's (The Mathworks, Natick, MA, USA) tri-cubic reconstruction kernels. The cubic kernel was chosen over linear interpolation to more closely resemble the ideal reconstruction kernel.

The original grid with (2.5mm) 3 voxel size was downsampled with factors k i = [1.5,2,2.5,3,4] to simulate the loss of image content by coarser imaging. After downsampling, the data was again upsampled by the respective inverse factor in order to be able to compare data on a common (the original) grid (see figure 1).

For the registration, baseline intensity variations throughout the imaging stack had to be removed. Accordingly, all images were equalized into the same dynamic range by normalizing the (magnitude) images to the median of the maxima found in all individual images. Subsequently, intensity values were capped to the estimated median-max. This was necessary because the T 1 -weighted imaging sequence yields the brightest signal for blood. Consequently, since the cardiac cycle has a similar time scale as the dynamic imaging, these signal maxima appeared modulated in vicinal vessels.

Motion estimation

Motion throughout the imaging session was estimated with an implementation [START_REF] Roujol | Real-time MR-thermometry and dosimetry for interventional guidance on abdominal organs[END_REF]) of non-rigid registration based on optical flow. As derived by [START_REF] Horn | Determining optical flow[END_REF], the optical flow's objective function is formulated as

Ω (I x v x + I y v y + I z v z + I t ) 2 + α 2 ∇v x 2 2 + ∇v y 2 2 + ∇v z 2 2 dxdydz (1)
The left side of the functional states that every spatio-temporal gray level variation ∇I and I t is attributed to motion v x v y v z T . In addition, the regularization term on the right side is introduced to penalize sudden spatial changes of motion, which was validated for organ deformations (Østergaard [START_REF] Noe | Acceleration and validation of optical flow based deformable registration for image-guided radiotherapy[END_REF]. Accordingly, the unknown DVF is then extracted by minimization of equation 1 with respect to

v x v y v z T .
In this work, the regularization parameter was set to α = 0.3. This value was chosen phenomenologically, in order to avoid instabilities in the motion estimation, that would lead to non-continuous estimations of the DVFs [START_REF] Roujol | Automatic nonrigid calibration of image registration for real time MR-guided HIFU ablations of mobile organs[END_REF].

Motion was estimated between the reference image (first image dynamic) and every consecutive imaging dynamic, .

✿✿✿✿✿✿✿

This

✿✿✿✿✿✿✿✿✿

yielded

✿ yielding 100 DVFs for the N = 100 dynamics. This process was repeated for every 4D stack of the respectively applied resampling factor with index i th .

Error assessment

As depicted in figure 2, the quality of the motion estimation for the i th dataset was assessed by taking the differences between the DVF from resampled data and the DVF I dyn (x,y,z,n) of the original data. However, the multidimensional, 4D vector fields

DVF i (x, n) =    v x v y v z    (x, n) (2)
require integral metrics for convenient quantification the differences. For clarity, the analyses are done solely in terms of absolute, component-wise differences (endpoint error [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF]))

∆DV F i (x, n) = DVF i (x, n) -DVF ref (x, n) 2 , (3)
with n denoting the temporal index of a 3D-dynamic.

Locally, two salient voxels were selected and analyzed over time

✿✿✿ for ✿✿✿✿ the ✿✿✿✿✿✿✿✿ longest ✿✿✿✿✿✿✿ dataset ✿✿✿✿✿✿✿✿✿✿✿✿ (volunteer ✿✿ 1 ✿✿ in ✿ table 1 ✿
). Furthermore, for generating global statistics, we make use of the root mean square (RMS) which is applied on the absolute component-wise differences between two DVFs

ǫ RM S,i (x) = 1 N N n=1 ∆DV F i (x, n).
(4)

Data selection

In order to limit the sample size, comparatively mobile voxels including continuous and spontaneous motion events were elected for the global statistics. Using linear regression, coefficients β(x) were calculated for every voxel over time as argmin

α(x),β(x) {α(x) + β(x)n = DV F ref (x, n) 2 } . ( 5
)
Voxels with β(x) > 0.02mm (dynamic) -1 were considered in the global statistics. In combination with the metric in eq. 5 this threshold extracts the elected dynamic voxels.

In addition to regression coefficients, the RMS of the DVFs could be a feasible quantifier of the target mobility in a voxel over time. However, for the chosen MR acquisition type, it did not prove feasible because of high false-positive rates. These were caused by the previously mentioned pulsation artifacts in the aorta and vena cava, randomly modulating the individual dynamics and thus the extracted DVF. Figure 3: Absolute DVF in a selected slice of volunteer 1 for increasing downsampling (a-f). Qualitative similarity is apparent until downsampling with factor 2.5.

Results

Figure 3 shows a selected slice of anatomy ✿ imaging dynamic for either of the 6 (resampled) datasets(-).

Over the selected slice, the DVFs extracted from the individually resampled data appear visually congruent. The stretched displacement at the edge of the liver, next to the right kidney, appears consistently throughout all datasets and holds its structure up to a two-fold downsampling before it disintegrates.

Large deformations in the region between stomach and duodenum are prominent in all scenarios and retain shape. However, starting from factor 2.5, a visible blur is introduced to the DVF distribution. At 4x-downsampling, formerly prominent spots significantly decrease in motion magnitude. In addition to the average error, higher moments of the error statistics increase equally, which leads to outliers ǫ RM S,i (x) > 2mm for downsampling with factors greater than 2.

In addition to the global statistics of many points (figure 5), local performance evaluation is shown in figure 6, illustrating the transient impact of downsampling. In the exemplarily selected point at the liver/kidney (figure 6a), a constant drift of about 2.5mm/17min is expressed over time. With coarser sampling starting from factor 3, the deformation appears greatly deteriorated with RMS errors proportional to magnitude. However, until factor 2.5, errors appear independent of the magnitude and stay lower than 1mm.

Slightly different in the intestinal sample (figure 6b), absolute errors start to exceed 1mm already at factor 2.5. Undersampling with factors 3, as before, results in a greatly deteriorated DVF with errors exceeding 2mm.

Discussion

Although removing image content by spatial downsampling, the used image registration algorithm was able to recover deformation fields to a high degree. In the synthetic experiments the most finely resolved data ((2.5mm) 3 ) served as ground truth; with this, congruent DVF were able to be extracted up to a downsampling factor of 2 ((5mm) 3 voxel size), with a local error of less than 1mm and global RMS deviations staying beneath 1.5mm. The boxplots illustrate the evolution of the RMS error within the elected voxels (β(x) > 0.02mm(dyn.) -1 ) for the various resampling factors and volunteers. In each cluster the RMS error statistics from volunteer 1, 2 and 3 is plotted from left to right. For completeness, a boxplot for the original, non-resampled case was added.

This observation contradicts the idea that only images of diagnostic quality are eligible for image registration purposes. Just the opposite, deformable image registration is much more data sufficient. Accordingly, imaging for the generation of displacement data is less demanding on image resolution than e.g. diagnostic imaging.

The fact that estimation of non-rigid motion requires spatial regularization is backing this hypothesis. Adding a spatial smoothness constraint (eq. 1, right side) enables stable motion estimation for regions with small image gradients [START_REF] Bruhn | Lucas/Kanade meets Horn/Schunck: Combining local and global optic flow methods[END_REF]. Additionally, to suppress the impact of outliers, literature suggests the application of low-pass filters similar to those used in this work for image preconditioning in order to obtain better motion estimations [START_REF] Barron | Performance of optical flow techniques[END_REF].

The simulated low-resolution MRI proposed here inherently implements both spatial (velocity) regularization and image low-pass filtering.

In real-time MRI, where image resolution can and must be traded off against imaging time, the proposed method could be used to significantly accelerate the acquisition. Using 2x spatial undersampling, the imaging time of about 3s for a (2.5mm) 3 -acquisition could ideally be sped up fourfold while retaining the quality of extracted DVF.

It is important to note that offline resampling using the cubic reconstruction kernels is solely an approximation of the effect of an actual low resolution MR acquisition, due to accompanying changes in contrast and point spread function (PSF). However, in first approximation, the used reconstruction kernels are methodologically sound. Combined with gated MRI, this technique enables fast 3D anatomy updates of abdominal volumes during the exhalation plateau. Used in an acquisition without respiratory trigger, i.e. free breathing, motion estimates would be less impacted by intra-scan artifacts (blurring) due to the shorter acquisition time.

Due to respiratory gating in the exemplary data, only medium-to long-term drift motion in the abdomen was monitored. Thus, observed deformations within the body are limited and spatially smooth to a first approximation. The image registration is constrained equally, having a spatial smoothness regularization.

Accordingly, similar results are expected with slowly changing pelvic sites. However, results may be different when observing steep gradients in DVF over time, such as when regarding abdominal breathing excursions. This would put stress on the differences between the algorithmic and the physiologic constraints. Due to the lack of (real-time 3D) data being eligible to be considered as ground-truth, this scenario has not been considered in this work. A potential data source could be retrospectively binned 4D MRI data.

Possibly, the accepted error could be evaluated per site or per subject with pre-treatment assessments: high-resolution data could be acquired and gradually downsampled until remarkable sites exceed an error threshold. These assessment could take place under worst-case treatment conditions, e.g. during rectal or bladder filling in pelvic cases, in order to test the error for an extreme scenario. The reward would be a highly accelerated MRI acquisition, without relying on the further introduction of signal models.

Conclusion

This investigation showed, that deformation estimation on coarsely resolved images yields congruent results when compared to more highly resolved images. The examined insensitivity of nonrigid image registration to spatial downsampling can be used as a facility to speed up MRI acquisitions, specifically for real-time therapy guidance. Careful optimization between acceptable registration error and necessary imaging speed will however be necessary, since this trade-off cannot be made in an intuitive way.

Figure 1 :

 1 Figure 1: Downsampled images resampled back into the original grid. Images exemplarily show the image information loss with higher downsampling factors

  ✿

Figure 2 :

 2 Figure 2: Evaluation workflow for motion estimation and error quantification: k i indicates the respective factors for the resampling operators (↓,↑), OF means the motion estimation using optical flow. Subsequent L 2 -norm ( • 2 ), ∆ and RMS operators work on a per-pixel basis. x,y,z are the spatial coordinates, n denotes the temporal index and i is the index for the individual downsampling factors.

  image stack. A transparent overlay is applied onto anatomy, representing absolute displacement fields estimated in the ✿✿✿ last ✿✿ (100 th )

Figure 4 :

 4 Figure 4: Global error quantification: distribution of the extracted regression coefficients β(x) over the anatomies of 3 volunteers

Figure 5 :

 5 Figure5: Global RMS error statistics: The boxplots illustrate the evolution of the RMS error within the elected voxels (β(x) > 0.02mm(dyn.) -1 ) for the various resampling factors and volunteers. In each cluster the RMS error statistics from volunteer 1, 2 and 3 is plotted from left to right. For completeness, a boxplot for the original, non-resampled case was added.

Figure 6 :

 6 Figure 6: Two manually selected points of interesting DVF evolution in volunteer 1 for multiple downsampling factors

Table 1 : ✿✿✿✿✿✿✿✿✿✿ Overview ✿✿✿ of ✿✿✿✿ the ✿✿✿✿✿✿✿✿✿ acquired ✿✿✿✿✿✿✿✿✿ amount ✿✿ of ✿✿✿✿✿✿ data ✿✿✿ for ✿✿✿✿✿ each ✿✿✿✿✿✿✿✿✿✿✿ volunteer. ✿✿✿✿ The ✿✿✿✿✿✿✿✿✿ amount ✿✿✿ of ✿✿✿✿✿ data ✿✿✿✿✿✿✿✿✿✿ acquired ✿✿✿✿ for ✿✿✿✿✿ each ✿✿✿✿✿✿✿✿✿✿✿ volunteer ✿✿ is ✿✿✿✿✿✿✿✿ stated ✿✿✿ in

 1 total ,

					✿✿✿✿✿✿✿✿✿✿ Volunteer ✿	#	✿✿ of ✿✿✿✿ 3D	✿✿✿✿✿✿✿✿✿✿ dynamics ✿	✿✿✿✿✿ total ✿✿✿✿✿ scan ✿✿✿✿✿✿✿✿✿✿ duration ✿
					✿ 1			100	✿✿✿✿✿✿✿ 17min
					✿ 2			✿✿✿ 150 ✿	✿✿✿✿✿✿✿ 15min
					✿ 3			✿✿ 80	✿	✿✿✿✿✿ 8min
								✿✿✿✿✿✿✿ From ✿✿✿✿ the
	✿✿✿✿✿✿✿✿✿✿ individual ✿✿✿✿✿✿ total ✿✿✿✿✿✿ scan ✿✿✿✿✿✿✿✿✿✿ duration ✿✿✿✿✿ and	✿✿✿✿ the ✿✿✿✿✿✿✿✿✿ amount ✿✿✿ of ✿ 3D dynamicswere measured over
	17min, with an approximate time per dynamic of 3s. The volunteer was	✿ , ✿✿ it	✿✿ is	✿✿✿✿✿✿✿✿✿ possible ✿✿ to
	✿✿✿✿✿ infer	✿✿✿✿✿✿✿ gating ✿✿✿✿✿✿✿✿✿✿✿ efficiencies ✿✿✿✿ for	✿✿✿✿✿ each ✿✿✿✿✿✿✿✿✿✿✿ volunteer:	✿✿✿✿✿✿✿ While	✿✿✿ for ✿✿✿✿ the ✿✿✿✿✿ first	✿✿✿✿✿✿✿✿✿✿✿ volunteer, ✿✿✿ an	✿✿✿✿✿✿✿✿✿✿✿✿✿ approximate
	✿✿✿✿✿✿✿✿ amount ✿✿✿ of	✿✿ 6	✿✿✿✿✿✿✿ images ✿✿✿✿✿✿✿ could	✿✿✿ be ✿✿✿✿✿✿✿✿✿✿ extracted ✿✿✿✿✿ per	✿✿✿✿✿✿✿✿ minute, ✿✿✿✿✿✿✿✿✿✿ scanning ✿✿✿✿ the ✿✿✿✿✿✿ other ✿✿✿✿✿ two	✿✿✿✿✿✿✿✿ yielded ✿✿✿ 10
	✿✿✿✿✿✿✿ images	✿✿✿✿ per ✿✿✿✿✿✿✿✿✿ minute, ✿✿✿ on ✿✿✿✿✿✿✿✿✿ avarage. ✿✿✿✿✿✿✿✿✿ During	✿✿✿✿ the	✿✿✿✿✿✿ scans ✿✿✿✿ the ✿✿✿✿✿✿✿✿✿✿✿✿ volunteers ✿✿✿✿✿ were ✿ breathing freely
	and was

✿ table 1 ✿ . ✿✿✿✿✿ were ✿ not commanded any breathing alteration.

Acknowledgements

The authors thank the ITEA (project 12026, SoRTS), the European Research Council (project ERC-2010-AdG-20100317, Sound Pharma) and Elekta AB (Stockholm, Sweden) for funding.


Markus Glitzner
email: m.glitzner@umcutrecht.nl

Sjoerd P M Crijns

Baudouin Denis De Senneville

Charis Kontaxis

Fieke M Prins

Jan J W Lagendijk

Bas W Raaymakers

On-line MR imaging for dose validation of abdominal radiotherapy

Keywords: respiratory motion compensation, dose reconstruction, MR guidance, adaptive radiotherapy, hypofractionated treatment, local deformations

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In an attempt for higher tumor control and better organ at risk (OAR) sparing, new motion management strategies are gaining interest.

Traditionally, geometric margin concepts such as planning target volume (PTV), internal target volume (ITV) (ICRU 2010) or mid-ventilation [START_REF] Wolthaus | Mid-ventilation CT scan construction from four-dimensional respirationcorrelated CT scans for radiotherapy planning of lung cancer patients[END_REF] approaches for uncertainty compensation are used. Additionally, various active motion compensation methods have been recently investigated . Amongst them are guided patient positioning [START_REF] Borst | Kilo-voltage conebeam computed tomography setup measurements for lung cancer patients; first clinical results and comparison with electronic portal-imaging device[END_REF]) as well as real-time beam gating [START_REF] Keall | The clinical implementation of respiratory-gated intensity-modulated radiotherapy[END_REF][START_REF] Shirato | Physical aspects of a real-time tumor-tracking system for gated radiotherapy[END_REF] and tracking techniques [START_REF] Keall | Motion adaptive x-ray therapy: a feasibility study[END_REF][START_REF] Keall | The first clinical implementation of electromagnetic transponder-guided MLC tracking[END_REF][START_REF] Crijns | Proof of concept of MRI-guided tracked radiation delivery: tracking one-dimensional motion[END_REF]. While the former are mostly X-ray based imaging techniques, the latter are controlled by feedback variables extracted from external or internal surrogates, which ideally move in phase with the target volume [START_REF] Ruan | Image-guided positioning and tracking[END_REF]. Obviously, this in-phase condition is non-trivial to maintain generally, as system lag and imaging uncertainties introduce position errors [START_REF] Keall | Geometric accuracy of a real-time target tracking system with dynamic multileaf collimator tracking system[END_REF][START_REF] Glitzner | On the suitability of Elektas Agility 160 MLC for tracked radiation delivery: closed-loop machine performance[END_REF]. Furthermore, using target surrogates for target tracking, motion of surrounding tissue in the beam's path is entirely neglected. Additionally, although several authors have previously investigated dose tracking using surrogate signals (e.g. [START_REF] Poulsen | A method of dose reconstruction for moving targets compatible with dynamic treatments[END_REF][START_REF] Ravkilde | Fast motion-including dose error reconstruction for VMAT with and without MLC tracking[END_REF])), the assumed (rigid) translations of the entire patient volume in phase with the surrogate signal are not physiologically sound.

The fusion of magnetic resonance imaging (MRI) diagnostics and radiotherapy treatment [START_REF] Lagendijk | The Magnetic Resonance Imaging-Linac System[END_REF][START_REF] Fallone | The Rotating Biplanar Linac-Magnetic Resonance Imaging System[END_REF][START_REF] Mutic | The ViewRay System: Magnetic Resonance-Guided and Controlled Radiotherapy Semin[END_REF][START_REF] Keall | The Australian Magnetic Resonance Imaging-Linac Program[END_REF] enables new options for motion management. With these combined modalitites, volumetric information is available on-line, i.e. during irradiation. This 3D information can serve for both target tracking (plan adaptation) and dose tracking (plan validation) [START_REF] Kupelian | Magnetic Resonance-Guided Adaptive Radiotherapy: A Solution to the Future Semin[END_REF]. However, although MRI offers unprecedented soft-tissue contrast, volumetric imaging times are in the order of seconds. In the abdomen, real-time MRI in 3D is thus technically challenging, due to the high mobility of the abdominal organs [START_REF] Langen | Organ motion and its management[END_REF][START_REF] Bussels | Respiration-induced movement of the upper abdominal organs: A pitfall for the threedimensional conformal radiation treatment of pancreatic cancer[END_REF] causing motion induced artifacts.

The aspiration of this work is to develop a technique to enable joint respiratory motion compensation and detailed 3D plan validation for the treatment of highly mobile abdominal organs using an magnetic resonance (MR)-linac modality. Therefore, the presented method covers two aspects, i.e.

• apply respiratory motion compensation using responsive, real-time radiation gating on an internal surrogate (diaphragm) and • recover the actual dose deposition by using 3D imaging over the active treatment time.

The unique combination of these fundamentals enables respiratory motion compensated, quality-assessed radiotherapy of abdominal sites. As a proof of concept, we evaluated two (virtual) kidney cases. This site was chosen because of its high mobility [START_REF] Schwartz | Kidney mobility during respiration[END_REF][START_REF] Stam | Kidney motion during free breathing and breath hold for MR-guided radiotherapy[END_REF] during respiration and its expected benefit from hypofractionated treatments [START_REF] De Meerleer | Radiotherapy for renal-cell carcinoma[END_REF]. Planned and true dose depositions were calculated and compared offline, integrating recorded treatment machine parameters, MR-acquired images and deformable image registration. Figure 1: Planning images (a and b) and intrafraction images (c). Kidneys, liver, myelum, pancreas and bowels were delineated in addition to the GTV (green) and PTV (red). In the depicted anatomy, a benign cyst was delineated as target volume in the kidney. Anatomic directions: head (H) and left (L).

Methods

In this paper, rapidly sampled MR-navigator echoes were used to sample the respiration and restrict both the treatment beam and the image acquisition to the quasi-static expiration phase of the breathing cycle (Section 2.1). During the gating interval, delivery took place and a fast 3D images were taken from the target region and its entire surrounding, i.e. all OAR in the beam's path. Due to the respiratory motion compensation, the images contained solely cycle-to-cycle and long-term variation, such as spontaneous and drift motion. Motion estimations (sections 2.2 and 2.3) extracted from the 3D imaging volumina are then used to recover anatomy variations in the course of treatment.

By continuous logging during delivery of the generated plan (sections 2.4 and 2.5) the individual anatomies could be correlated with the machine state at every instance. Subsequently, applying Monte-Carlo (MC) dose calculations (section 2.7) on the pseudo-computed tomographys (CTs) of all individual anatomies (section 2.6), the actual delivered dose was reconstructed. In this way, an iterative dose accumulation was performed over the entire treatment time.

MRI acquisition

All imaging experiments were performed on a 1.5T MRI machine (Ingenia, Philips Medical Systems, Best, The Netherlands) using a 28-channel torso array.

Two different MR pulse sequences were developed. Both of them employed navigator echoes [START_REF] Wang | Navigator-echo-based realtime respiratory gating and triggering for reduction of respiration effects in three-dimensional coronary MR angiography[END_REF] for respiratory gating in the exhale phase of the respiratory cycle. The navigator beam was placed on the right liver dome in order to maximize sensitivity to breathing excursions in caudo-cranial (CC) direction.

The images from the first pulse sequence were used for generating the delineation reference (figure 1a and1b). The sequence consisted of navigator-gated, balanced turbo field echo (TFE) shots which were triggered in a narrow gating window of 2mm (exhalation) in order to maximally suppress breathing artifacts. The shot length was approximately 250ms including an α/2-instance and 10 additional startup cycles to equilibrate the MR signal.

For fast 3D imaging during treatment a second sequence was developed. It consisted of a navigator echo (5mm gating window at exhalation) followed by a binomial, water-selective pulse. Subsequently, a multi-shot echo planar imaging (EPI) readout was driven to increase sampling efficiency due to the long navigator and water-selective preparation. The long water-selective pulse was chosen to enhance the performance of the image registration algorithm used (Section 2.2) which greatly benefited from suppressed fat signals.

The EPI readout is prone to geometric distortions due to susceptibility variations [START_REF] Haacke | Magnetic Resonance Imaging: Physical Principles and Sequence Design[END_REF]. Thus, the EPI train length was limited to 21 k-lines per shot to lower image acquisition times to around 3s. Additionally, moderate sensitivity encoding (SENSE) was employed with undersampling factors of P = 1.5 and S = 1.5 in the anterior-posterior (AP) and CC phase encoding direction, respectively.

The sequence parameters were T R /T E /α = 47ms/9.8ms/25 • . Typical gradient characteristics reported by the system were 192T /m/s for the slew rates, 12mT /m EPI-blip strength and 22mT /m readout strength.

As volume of interest (VOI), a 10cm thick transverse slab was selected (figure 1c) to cover the lower liver/upper right kidney and thus the treatment volume and its surroundings with (2.5mm) 3 isotropic resolution.

Motion estimation

The isotropic 3D data MRI acquired while treatment is inherently suited for non-rigid image registration. In this work an implementation of optical flow-based registration [START_REF] Roujol | Real-time MR-thermometry and dosimetry for interventional guidance on abdominal organs[END_REF]) was employed. The optical flow's objective function by [START_REF] Horn | Determining optical flow[END_REF] reads,

Ω (I x v x + I y v y + I z v z + I t ) 2 + α 2 ∇v x 2 2 + ∇v y 2 2 + ∇v z 2 2 dxdydz (1)
with I t , I x , I y and I z being the temporal and spatial derivatives. Next to the motion estimation using image gradients (left side of (1)) a data regularization term with a single parameter α penalizes sudden spatial changes of the estimated motion (∇v x , ∇v y , ∇v z ,). This regularization has been introduced for therapy guidance on mobile organs [START_REF] De Senneville | Real-time adaptive methods for treatment of mobile organs by MRI-controlled high-intensity focused ultrasound[END_REF]) and validated for organ deformations (Østergaard [START_REF] Noe | Acceleration and validation of optical flow based deformable registration for image-guided radiotherapy[END_REF]. Next to the single parameter, a further advantage is the performance of the optical flow algorithm when implemented on graphics processing units (GPUs) [START_REF] Roujol | Real-time MR-thermometry and dosimetry for interventional guidance on abdominal organs[END_REF]. In this work the regularization parameter α, was set to a conservative α = 0.3, in order to avoid non-continuous motion estimates [START_REF] Roujol | Automatic nonrigid calibration of image registration for real time MR-guided HIFU ablations of mobile organs[END_REF]. Motion estimation of the continuously acquired volumes (144px × 144px × 36px) took approximately 230ms using a compute unified device architecture (CUDA) implementation on an Geforce GTX280 (NVIDIA, Santa Clara, CA, USA). The host computer was equipped with 2 central processing units (CPUs) (Penryn quad-core, Intel, Santa Clara, CA, USA) with 8GB DDR3 SDRAM memory.

Gridding strategies

Generally, motion is estimated between two dynamics: One is the reference (I ref ) holding the initial anatomy. The second (I n ) is a version of the reference, which underwent physiologic deformation.

When estimating motion, generally, a pull-back deformation vector field (DVF)

V i (x) is calculated for each dynamic image I i (x), so that reference V n (V n ) -1 n=1 n=N (a) DVF estimation pCT (V n ) -1 n=1 n=N (b) pCT generation V n Σ n=1 n=N (c) accumulation
I ref (x) ≃ I n (x + V n (x)) .
(2)

This projects the content of the n th dynamic's grid to the reference grid. In this work, this transformation is consequently used to warp dose deposited on an intermediate grid back onto the reference. Is it accumulated on the reference grid (figure 2c).

Contrary, the push-forward DVF is transforming content from the reference grid onto the dynamic grid, i.e.

I n (x) ≃ I ref x + V -1 n (x) . (3)
Herein, as illustrated in figure 2b, this transformation is employed to construct dynamic CT data corresponding to the anatomical changes measured by the continuous 3D MRI as proposed by [START_REF] Boye | Mapping motion from 4D-MRI to 3D-CT for use in 4D dose calculations: a technical feasibility study[END_REF].

Figure 3 depicts the difference between push-forward and pull-back DVF. The latter can be estimated by interchanging reference and dynamic image in the motion estimation algorithm. Ideally, V n and V -1 n are inversely consistent [START_REF] Yang | A fast inverse consistent deformable image registration method based on symmetric optical flow computation[END_REF]) such that

V n (x + V -1 n (x)) + V -1 n (x) = 0. ( 4
)
In practice however, inverse consistency is generally not achieved by motion estimation algorithms. Thus, in this paper, inverse consistency was approximated by an iterative solution previously described in [START_REF] Heinrich | MIND: Modality independent neighbourhood descriptor for multi-modal deformable registration[END_REF] after calculating both push and pull DVF. Table 1: Planning constraints and parameterization for the created kidney case.

Plan generation

For this proof of concept, renal cell carcinoma (RCC) cases were generated by projecting fictional lesions into two volunteers' right kidneys. Additionally, both kidneys, myelum, liver, stomach, pancreas and bowels were identified as OAR and contoured as seen in figure 1a.

Based on constraints previously published by Stam, van Vulpen, Barendrecht, Zonnenberg, Crijns, Lagendijk & Raaymakers (2013a), a plan was generated using the parameters in table 1a. For this study, the planned isocenter was set to the center of the tumor. The planning goal was to deliver 25Gy to the target. A PTV-GTV margin of 2mm was used to obtain the 80% isodose line at the border of the GTV (D 99,GT V = 20Gy). Fluence was optimized by in-house developed software using an iterative approach by [START_REF] Ziegenhein | Performance-optimized clinical IMRT planning on modern CPUs[END_REF] using pre-calculated beamlets for 15 equally space beam angles (see table 1b). The MC-based beamlet generation [START_REF] Hissoiny | Fast dose calculation in magnetic fields with GPUMCD[END_REF] took the MR-linac magnetic field B 0 = 1.5T into account. Finally, the plans were sequenced into deliverable stepand-shoot intensity modulated radiotherapy (IMRT) segments using the adaptive sequencing (ASEQ) method [START_REF] Kontaxis | Towards adaptive IMRT sequencing for the[END_REF].

Plan delivery

Subsequently, the generated plan was delivered on a clinical linac emulator (Elekta AB, Stockholm, Sweden). The emulator featured an Agility 160 multileaf collimator (MLC) (Elekta AB, Stockholm, Sweden) and a realistic model of the gantry and beam generation. During the entire delivery the machine parameters were logged with a sampling interval of 40ms. This step was required to synchronise the exact machine status (leaf position, diaphragm position, gantry angle, dose rate) to the respectively valid anatomy.

Pseudo-CT definition and generation

Generally, for the planning stage, information about the electron density is needed for the fluence optimization. However, since CT data acquisition of the (healthy) volunteers was not feasible, pCTs were synthesized. Based on a two-component model (Kerkhof et al. 2010, Stam, van Vulpen, Barendrecht, Zonnenberg, Crijns, Lagendijk & Raaymakers 2013a), the body outline was extracted from the reference MR image. Subsequently the body was filled with the Hounsfield units (HUs) of water HU water = 0, while the surrounding air was set to the HU of air HU air = -2900. This yielded the two-component pCT.

Subsequently, the push-forward DVF V -1 extracted from the dynamic data (figure 2a) was used to transform the planning-pCT to the new dynamic grid (figure 2b). This yielded a pCT for every dynamic, pCT n .

Dose reconstruction

In this work, the navigator echo signal was used to trigger both imaging and the simulated treatment beam. For every pCT instance, an integration time ∆T n is calculated by summing the valid gating intervals as depicted in Figure 4. Together with the machine parameters extraction described in section section 2.5, a relation between momentaneous anatomy and machine state can be built up. Every acquired 3D volume is thus valid for one or more individual segment(s).

Subsequently, using a MC based dose calculation [START_REF] Hissoiny | GPUMCD: A new GPU-oriented Monte Carlo dose calculation platform[END_REF], dose is reconstructed on the respectively valid pCT n using the individual segment from the logged machine data [START_REF] Luo | Monte Carlo based IMRT dose verification using MLC log files and R/V outputs[END_REF]. This produced dose maps D n for all acquired anatomies, which were eventually warped back to the original grid using V n (figure 2c).

For MC calculations, a computer with two 12-core (Intel Xeon E5-2695, Intel, Santa Clara, CA, USA) CPUs and a Tesla K20c (NVIDIA, Santa Clara, CA, USA) GPU was used. One segment took around 15s to calculate with a 5% variance.

Plan validation

In order to validate the treatment plan, nominal (planned) doses were compared against the reconstructed doses. For error quantification, dose-volume histogram (DVH) points and integral parameters (D max , D mean , D min ) were calculated for the target volumes and the surrounding organs at risk. Furthermore, planning constraints were compared between planned and actual delivery. To pinpoint local differences, dose difference maps were calculated and the respective maximal, mean and minimal deviation was extracted. Figure 5 depicts representative transverse slices of the anatomic reference images for both cases. The kidney lesions differ in location: case 1 is located in a more internal, inferior position as compared to case 2. The visual comparison between the planned dose on the rigid anatomy and the reconstructed dose based on variable anatomy shows no obvious deviation in the selected slice. The global DVH-points shown in table 2 show similar features: all dose constraints defined in table 1a were satisfied both with the static planning case and after validation using the variable anatomy. Only kidney-D 25 of case 2 is close to its maximum dose of 5Gy. However, neither in the nominal nor in the reconstructed dose, the constraint is significantly violated and is therefore considered as plan compliant.

Results

A

For case 1 local differences in table 2a shows little variation when incorporating moving anatomy. This is consistent with the global DVH point for this case, showing only little differences between static and moving anatomy.

Regarding the parameters for the second case (table 2b) reveals discrepancies between the DVH points and the local differences. While the differences in PTV are only around 1%, local underdosage exceed 2Gy with overdosages up to 1.5Gy in the PTV.

As the most vicinal OAR, the statistics for the liver show a similar pattern. Although, DVH points do not suggest changes in the actually deposed anatomy, locally variations are significant, ranging from -1.3 to 1.9Gy.

These differences in the liver-kidney region are equally prominent in the difference images of figure 6b. As primary source for the dose alteration, the beam at 216 • (posterior-anterior ) induces high difference gradients in the kidney-liver transition region. Both underdosage and overdosage in this region is caused by anatomy variations in this beam angle. These variations are especially visible in the transverse and coronal dose profiles in figure 6b and are correlated with the internal motion between kidney and liver. At the exiting portal of the same beam, the exit dose changed due to a variation of the diaphragm. Figure 6a as a comparison shows the same exit dose effects for the posterior entering beams. However, effects in the vicinity of the GTV are minor.

Discussion

Radiotherapy of mobile sites requires treatment validation. This work pursues a technique for on-line dose tracking for the treatment of abdominal sites using the on-line image capabilities of an MR-linac. We successfully devised and implemented a treatment pipeline capable of reconstructing the dose deposition in kidney cases using on-line MR imaging during treatment. By use of fast MR-navigator echoes, the abdomen was decoupled from breathing excursions. Entering the gating window triggered both treatment beam and 3D image acquisition, the image covering an entire transverse slab around the target volume.

The continuous 3D acquisitions can be utilized manifold: Firstly, deposited dose can be calculated integrally for an entire fraction. This enables the option to adapt subsequent fractions according to the intermediate deposition. Also, particularly erroneous parts of delivery causing high dose differences can be traced back in retrospect. Delivery elements which turn out to systematically induce larger errors in the plan validations can be identified and avoided for the sensitive population.

Secondly, when performing real-time accumulation after each imaging volume, local dosimetric or geometric errors at OAR can be used to trigger exception handling, based on predefined levels.

The two cases processed in this work suggest that inter-patient variations are significant. While for case 1 no considerable dose differences were found locally, the second case revealed considerable local differences. Most remarkably, the integral DVH points did not significantly change at any time. This suggests that integral measures (γ, DVH points) alone are insufficient to characterize the validity of the plan outcome in a dosimetric sense.

Evidently, using on-line 3D MR imaging for the dose reconstruction resolves the necessity to model anatomic motion as a bulk translation according to a surrogate signal, as described by [START_REF] Poulsen | A method of dose reconstruction for moving targets compatible with dynamic treatments[END_REF].

Compared to electronic portal imaging device (EPID) dosimetry as described for example by [START_REF] Lin | 4D patient dose reconstruction using online measured EPID cine images for lung SBRT treatment validation[END_REF], the presented method resolves the entire beam path, rather than projecting back over an integral pathway. Additionally, projective dosimetry relies on previously (retrospectively) acquired patient anatomies. Using EPIDs in the static B-field of an MR-system, the dosimetry would heavily rely on MCgenerated models of the beam's trajectory, in order to account for electron return effects (EREs) [START_REF] Raaymakers | Integrating a MRI scanner with a 6 MV radiotherapy accelerator: dose increase at tissue-air interfaces in a lateral magnetic field due to returning electrons[END_REF]) at tissue-air boundaries.

In contrast, with the proposed pipeline, anatomy information is available on-line. Using the gated approach, anatomical drifts and other cycle-to-cycle variations, e.g. shown by Seregni et al. (2015), do not corrupt the anatomical accuracy. Instead of inferred data from a retrospective 4D model as discussed by [START_REF] Sawant | Investigating the feasibility of rapid MRI for image-guided motion management in lung cancer radiotherapy[END_REF], true 3D data is sampled. This yields accurate data representing the momentary state of anatomy.

For a potential clinical deployment, the influence of target motion within the gating window has to be taken into account as shown by Stam, van Vulpen, Barendrecht, Zonnenberg, Crijns, Lagendijk & Raaymakers (2013b). This can be accomplished by trading off gating window width against a residual motion margin.

Since the EPI readout in the imaging sequence is prone to geometric distortions, investigations of distortions correction methods will be necessary. The latter were already addressed by, e.g., [START_REF] Dragonu | Real-time geometric distortion correction for interventional imaging with Echo-planar imaging (EPI)[END_REF] for real-time high intensity focused ultrasound (HIFU) interventions. Moreover, to overcome image distortions, highbandwidth fast field echo (FFE) sequences will be considered, despite their contrast and artifact behavior, which have an adverse effect on the image registration.

The gating efficiency was estimated to around 50% for both cases, which would translate into an approximately two-fold treatment time. To address the problem of gating window drifts, we are currently evaluating gating level adaptation schemes. Additionally, ways to validate the dose reconstruction workflow itself are currently under investigation at our institution. Furthermore, ongoing efforts are translating the current off-line implementation into an on-line pipeline, deployable on the MRlinac.

Conclusion

We successfully implemented and demonstrated a workflow to perform a respiratory motion compensated, fully reconstructed treatment for abdominal site using the online imaging capabilities of an MR-linac. This proof of concept is a step towards the clinical implementation of dose validation using such a device.

  Reference with delineations (192px × 192px × 100px, (2mm) 3 ) (b) Reference with treatment VOI (192px × 192px × 100px, (2mm) 3 ) (c) VOI during dynamic imaging (144px × 144px × 36px, (2.5mm) 3 )

Figure 2 :

 2 Figure 2: DVF estimations (a) and regridding for pCT-generation (b) and dose accumulation (c).

  pull-back DVF V push-forward DVF V -1

Figure 3 :

 3 Figure 3: warping between grids; The pull-back DVF provides the mapping back to initial grid of the circle. The inverse, V -1 , serves to reconstruct values from the circle's grid onto the ellipse's grid. Deformed indexing grids are shown in blue. Arrows indicate the direction of registration. Site Parameter Value GTV D 99 20Gy Kidney D 25 ≤ 5Gy D mean ≤ 8Gy Liver D 30 ≤ 2.5Gy Bowels D 5 ≤ 20Gy D max ≤ 21Gy Stomach D 4 ≤ 22.5Gy Myelum D max ≤ 13Gy (a) Planning constraints

Figure 4 :

 4 Figure4: Navigator data with acceptance window: between the gating levels, imaging and beam generation becomes active. While the MR volumes are sampled, the plan segments are gradually delivered according to the treatment plan and the gating intervals. The high gating efficiency (≈ 53%) is clearly recognizable.

Figure 5 :

 5 Figure 5: Cases 1 and 2 in figures (a) and (b), from left to right: selected slices (top row: transversal, bottom row: coronal) in delineated anatomy, the planned dose distribution and the accumulated dose incorporating the moving anatomies. Anatomic directions: anterior (A), head (H) and left (L).

Figure 6 :

 6 Figure 6: Visual dose and image comparison of both cases (a and b). Left: anatomic orientation, middle: dose differences between reconstructed delivery and planning (D recon -D plan ), right: difference between final and reference anatomy (I end -I ref ), indicating the anatomy displacements over the treatment time transversally (top row) and coronally (bottom row).

Table 2 :

 2 DVH points and difference statistics for treatment validation of both cases. For each VOI, DVH-points (D xx ), spatial integral parameters of absolute doses (D max , D mean , D min ) and voxel-wise differences (∆D max , ∆D min , ∆D mean ) are shown. Planning constraints are marked with an asterisk (*).

Acknowledgements

The authors thank the ITEA (project 12026, SoRTS), the European Research Council (project ERC-2010-AdG-20100317, Sound Pharma) and Elekta AB (Stockholm, Sweden) for funding.