M Glitzner
email: m.glitzner@umcutrecht.nl

S P M Crijns

B Denis De Senneville

J J W Lagendijk

B W Raaymakers

On the Suitability of Elekta's Agility 160 MLC for Tracked Radiation Delivery: Closed-loop Machine Performance

Keywords: MLC, tracking, system response, delay analysis, band-limited, linear-phase, physiologic motion

For motion adaptive radiotherapy, dynamic multileaf collimator (dMLC) tracking can be employed to reduce treatment margins by steering the beam according to the organ motion.

Until now, the Elekta Agility 160 MLC has hitherto not been evaluated for its tracking suitability under physiologic conditions. Both dosimetric performance and latency are key figures and need to be assessed generically, independent of the used motion sensor.

In this paper, we propose to use harmonic functions directly fed to the MLC to determine its latency during continuous motion. Furthermore, a control variable is extracted from a camera system and fed to the MLC. Using this setup, film dosimetry and subsequent γ statistics are performed, evaluating the response when tracking magnetic resonance imaging (MRI)-based physiologic motion in a closed-loop.

The delay attributed to the MLC itself was shown to be a minor contributor to the overall feedback chain as compared to the impact of imaging components such as MRI sequences. Delay showed a linear phase behaviour of the MLC employed in continuously dynamic applications, which enables a general MLC-characterization. Using the exemplary feedback chain, dosimetry showed a vast increase in pass rate employing γ statistics.

In this early stage, the tracking performance of the Agility using the test bench yielded promising results, making the technique eligible for the translation to tracking using clinical imaging modalities.

Introduction

In an attempt to further limit radiation exposure of healthy tissue, real-time tracking has been evolved recently. Potentially, using tracking, the planning target volume (PTV) can be decreased significantly. Treatment margins added to account for set-up and motioninduced misalignments can be compensated using image guided radiotherapy (IGRT), while dose-painting the actual tumor volumes retaining high conformity with intensity modulated techniques (Giraud et al. 2006). A review of image guided radiotherapy can be found in (Ruan et al. 2011).

For multileaf collimator (MLC) tracking, several concepts have been demonstrated for main manufacturers of linac hardware (Tacke et al. 2010, Falk et al. 2012). Multiple position extraction techniques are utilized to assess tumor motion indirectly via surrogates, or directly via fiducials.

In order to obtain geometric organ position, mostly surrogate based positioning has been in use, amongst them optical positioning systems (Sawant et al. 2009) and electromagnetic transponders (Smith et al. 2009) (Real-time Position Management (RPM) system and Calypso, Varian Medical Systems Inc., Palo Alto, CA, USA).

Recently, the first clinical application utilized the latter modality on a prostate case (Keall et al. 2014). Moreover, tracking signals have been extracted from kV and MV imaging data (Berbeco et al. 2005, Cerviño et al. 2009) and recent studies also investigated dMLC tracking with the use of real-time MRI-signal (Crijns, Raaymakers & Lagendijk 2012, Yun et al. 2013).

To this day, the performance of real-time MLC tracking with Agility 160 MLCs (Elekta AB, Stockholm, Sweden) has not yet been explored. However, Davies et al. (2013) reported on virtually tracked volumetric intensity modulated arc therapy (VMAT) using initial synchronization between the delivery and a motion platform.

The aspiration of this paper is to assess the system performance of the MLC in a generic sense. For the application in an integrated tracking modality, knowledge of the time delay of an MLC is fundamental, as it is the most sensitive parameter determining the geometric error in motion compensation (Poulsen et al. 2010, Ruan 2010), parameterization of prediction filters and retrospective dose accumulation. The relative contribution of the MLC is needed for evaluating its feasibility in a feedback chain with clinical components such as in an MR-Linac system developed in our institution (Raaymakers et al. 2009, Lagendijk et al. 2014). Due to the bandlimited character of physiologic motion, the time shift between aperture application and mechanic execution is determined in a physiologically reasonable continuous motion regime. To legitimate a general definition of the delay, the MLC is tested for linear phase behaviour in this range of harmonics.

Additionally, an optical feedback chain consisting of fast imaging, image processing and 1D aperture adaptation is prototyped which serves as a test bench for the extraction of timing characteristics of the MLC and for actual dosimetry. It can easily be adapted to other control variables, potentially sourced from other imaging techniques such as from magnetic resonance imaging (MRI), electromagnetic transponders or optical positioning surrogates. The final endpoint being radiation depositions in clinical applications, the dosimetric impact of the actual feedback loop is measured in response to a physiologic motion scenario. Tracked and untracked dose maps are compared with integral quantifiers, namely tight-radius γ-analysis and relative dose histograms, in order to rate the tracking performance without regarding for uncertainties such as from planning or from positioning.

Methods

The methods are divided into a general part comprising hardware and evaluation (sections 2.1 to 2.4) together with the description and measurements of the elements involved in open-and closed-loop MLC control (sections 2.5 to 2.7).

MLC characteristics

For all experiments, Elekta's Agility 160 MLC as part of a clinical Axesse (Elekta AB, Stockholm, Sweden) linear accelerator is used. The Agility features 80 leaf pairs and dynamic diaphragms to limit the aperture perpendicular to leaf travel direction. For reasons of data presentation and evaluation, an aperture of 20 × 40mm 2 was applied, by separating all leaves by 40mm. The arising gap is then limited to 20mm using the diaphragm. The leaf motion is in caudo-cranial (CC) direction, parallel the longer side of the rectangular aperture. The head is thus rotated to 90 • . The aperture width of 40mm is chosen to yield a global dose plateau for each measurement, i.e. a point of permanent irradiation. This is needed to obtain a reference measurement for the relative dosimetry, detailled in section 2.3.

Tracking components

The feedback loop comprised an optical path for object detection and a control loop which shifts the aperture in response to object motion. The optical path located the position of a binary image (target) mounted on a linearly (1D, along leaf motion direction) moving object. The target image contained a contrast-rich structure with sharp gradients (see figure 1b) in all directions in order to make it locatable by the image processing explained below. The camera was positioned in such a way, that the target image was moving in a fixed imaging plane perpendicular to the camera's optical axis. Secondly, the control loop polled the extracted object's position continuously and applied it in form of a rigid aperture translation (again 1D in leaf motion direction) to the MLC system.

The in-house developed tracking software served two purposes, depicted in figure 1: image acquisition and processing and commanding the MLC. In order to obtain the position of the tracked object, a camera (C930e, Logitech, Apples, Switzerland) presents a videostream (640px × 480px × 30Hz) containing the marked object and passes the stream to the image acquisition module. The camera was accessed by the computer vision library OpenCV (Bradski 2000). In order to obtain a reliable position extraction of the marked object by the subsequently described tracking algorithms, auto-focus features are disabled on camera-side.

In feedback experiments, the captured images are passed to a processing module without further conditioning.

The processing module extracts the pixel location of the marker based on a fast in-house minimum output sum of squared error (MOSSE) filter implementation (Bolme et al. 2010[START_REF] Crijns | Proc. Int. Soc. Magn. Reson. Med.'[END_REF], tuned for working in a 64 × 64px local region. Initially, before each experiment, geometric reference is established by a homing procedure in combination with the motion stage described below and a mean-free cross-correlation filtering. The latter extracts the position of an object marked with a known pattern from the captured images. Because of the relative slowness of this tracking method compared to the MOSSE filter, the correlation is only applied during homing. This algorithm transistion is necessary due to the local nature of the MOSSE filter; the detection region needs to be determined by the global cross-correlation algorithm first. In order to compare the processing performance, the runtimes of both algorithms were measured. Both filters operated on pixel pitch of 0.3mm/px. The change in object position extracted by the image processing module were used as the control variable for the 1D closed-loop aperture translation and the leaf positions of the translated aperture were applied to the Agility MLC control system via Ethernet with a cycle time of 40ms. Concurrently, the leaf positions provided by the MLC control system were continuously recorded.

All software was run via a Debian operating system on a computer equipped with 8GB of memory and a 3.4GHz quad-core CPU (Core i5-2400, Intel, Santa Clara, CA, USA).

Dosimetry

The dosimetric platform is sketched in figure 1 showing a film sandwiched between two layers of water-equivalent solid build-up material of 2cm thickness and a square surface of 10×10cm 2 each. Relative dosimetry is performed using the radiosensitive GafChromic EBT3 film (Ashland Inc., Wayne, NJ, USA), with dimensions of about 10 × 4cm 2 . The longer side of the film is parallel to leaf travel direction and thus to the direction of motion.

Generally, the film set-up moved in a plane of 1m SAD at isocenter. The object holding the film is referenced relative to the isocenter. In other words, an object in initial position (x = 0 in figure 2a) is under an isocentric aperture. With each irradiation, a dose of D = 1000M U is delivered using a constant dose rate of Ḋ = 400M U/min at 6M V .

For planar (2D) dosimetry, the irradiated film is scanned using an Epson 10000XL flatbed scanner (Seiko Epson Corporation, Suwa, Japan). An automatic alignment method centers and rectifies the dose images using rigid transformations.

For 2D evaluation, γ-statistics (Low et al. 1998) and dose difference test were used. The 2D dose distribution are transformed to relative dose maps ranging between noise floor and dose plateau. To establish a stringent quantifier for the synthetic experiments, a tight γ-pass radius of (1%/1mm) was set. The dose-window to be evaluated is between 10% and 100% to reject low-dose regions, which we consider unimportant for rating the tracking performance.

Phantom motion

Motion is induced using a home-built 1D translating stage. The linear positioning system (LPS) consists of an MDrive 17 (Schneider Electric Motion, Marlborough, CT, USA) driving a stage via a threaded rod (figure 1, green elements). The controller and drivers of the MDrive 17 are integrated and controllable via RS422. The motion stage carried the dosimetry set-up (cf. section 2.3) along with the optical target, detectable by the tracking algorithm (see figure 1b). In-house developed software was used to feed the stage with positioning data. The positioning data in the experiments were of two kind: synthetically generated sinusoidal curves and motion extracted from a temporally resolved liver position over multiple breathing cycles (see figure 2).

Sinusoids are applied to the motion stage for determining the overall closed loop response delay. Staying within a pseudo-continuous domain for the tracking system (the control variable does not strongly feature discrete step-behaviour), a parameterization of (A, T ) = (15mm, 6s) was chosen, where A is the amplitude and T the duration of a sinusoid. Physiologic motion patterns were extracted from CC components of deformation vector fields (DVFs) in the liver, extracted from coronally acquired MRI with 10Hz imaging frequency. The dosimetric response was evaluated based on these motion patterns in order to yield physiologically characteristic results. In order to present a continuous coordinate update to the motion controller, the data is upsampled by a factor of 2 to 20Hz using Matlab's (The Mathworks, Natick, MA, USA) standard reconstruction kernel (linear-phase Kaiser window filter). Consecutively, due to the relatively high noise level in the DVF, a 2.5Hz-Butterworth low-pass (n = 10) is applied to condition the raw motion signal before application to the motion stage. The noise level is caused by the echo planar imaging (EPI) nature of the acquisition and the consecutive DVF extraction by an algorithm using an optical flow signal model (Horn & Schunck 1981). The cutoff frequency of the filter is chosen to leave physiologic components of the spectrum untouched. Thus, it solely filters higher-frequency noise induced by the imaging and motion estimation. The strong drop-off of the motion spectra over approximately 0.3Hz in figure 2b is an indicator for the frequency cut-off of the sampling field in the open-loop benchmark (cf. section 2.7).

Lag determination

The retardation of the entire system is broken down to its elements illustrated in figure 3. Two major experimental pathways have been followed throughout the experiment. Experiments regarding the MLC as an isolated element used the open-loop configuration, where the aperture is directly applied to the MLC (see figure 3 lower solid path). In contrast, performing closed-loop experiments (figure 3 upper solid path), the MLC was steered by the feedback control variable extracted by the tracking software.

From a control input on the left side of the schematic in figure 3, the signal is sequentially passed through the respective processing steps which take up their individual times ∆T 1,2,3,4 . These elementary delays are regarded in detail subsequently in sections 2.6 and 2.7.

Closed tracking loop

For the delay determination of the overall tracking system, the mean time elapsed between detection of the target and the aperture reaching that target has to be evaluated. This chain consists of several abstract elements which are specified in the following.

Regarding again figure 3, object motion has to be captured by the camera hardware. After the readout, data is buffered and streamed via USB and the camera driver module until it is passed to the OpenCV application programming interface (API). The arising time delay is called ∆T 1 . The received image is then processed subsequently in order to extract position information with the filtering algorithms presented in section 2.2. The runtimes of these filters, from the instance of grabbing the image until presenting the extracted position information to the MLC control system, are denoted ∆T 2 . After exposing the current target position to the software interface, the MLC system needs time to approach the desired position. This lag is caused by the runtime of the MLC control system and the electromechanical properties of the MLC and is called ∆T 3 .

Finally, for logging purposes, data is being received from the MLC side. Like the object position, the MLC position logging is delayed, taking ∆T 4 to arrive at the control computer.

For measuring the delay of the closed tracking loop, harmonics of the form

x (t) = A sin 2π T (t + δ) (1)
are applied to the motion stage (cf. section 2.4) and tracked by the feedback loop. Logging both position updates and actual leaf positions yields ∆T loop = δ obj.pos -δ M LC,logged .

(2)

Since the position updates arrive and leaf positions arrive delayed by ∆T 1 + ∆T 2 and ∆T 4 , respectively, the loop latency is to be corrected with ∆T ′ loop,corr = ∆T loop + (∆T 1 + ∆T 2 -∆T 4 ) .

(3)

Identification of lag elements

Subsequently, the measurement of the system elements contributing to the overall delay of the feedback system are discussed.

Optical line delay ∆T 1 is the integrated lag value between a physical change of lighting and its detection in software. Offering transistor-transistor logic (TTL) outputs, a green light emitting diode (LED) is connected to the parallel DB-25 interface. In two separate threads, the light is switched on irregularly and measured by means of (binary) light detection using the same image acquisition software as in the tracking application. Timestamping both incidents yields a mean delay between incident and detection (Schellen et al. 2013). The parallel port features a sufficiently fast reponse time in the µs-range (Koolwal 2009), sufficing the immediate switching requirement of this global time referencing step. The LED's electric dynamics is considered to be delay-free.

Acquisition and processing delay ∆T 2 is the time it takes to grab an image from the OpenCV API and process and forward it to the MLC control. As the value is local runtime, the time ∆T 2 is directly logged in each pass.

MLC delay or open-loop lag is the latency of the MLC system only. It is determined by commanding the aperture by a single-threaded program running at a fixed cycle time of 40ms which corresponds to the real-time frame of the Agility MLC control system. Commands are sent to the MLC in each cycle. At the time of sending, the commanded position is recorded, as well as the actual leaf position retrieved by the logging stream introduced in section 2.2.

To quantify the phase of the leaf positioning, an amplitude-period (A, T ) parameter field of harmonics is sampled by applying synthesized harmonics directly to the MLC. The sines are in the form of 1, with parameter A being in the domain {1, 1.5, 2, 2.5, 3}cm and T in {3, 4, 5, ..., 8}s. Yielding pairs of sampled sine curves, the commanded and the observed, both curves are consequently fitted to a model as in 1 using matlab's Curve Fitting Toolbox. To obtain the dynamic MLC lag, the temporal displacement δ for both fits are related; as the observed state is delayed for the logging/subscription lag ∆T 4 , the measured phase difference is corrected for the subscription delay time T 4

∆T 3 = δ M LC,commanded -δ M LC,logged -∆T 4 . ( 4
)
Logging delay ∆T 4 is indirectly measured by capturing the aperture light field (ALF) in an open-loop experiment. The ALF position is extracted using the tracking method descibed in 2.2. ∆T 4 is then determined by fitting the ALF-mediated aperture position against the subscribed positions. The influence of ∆T 1 and ∆T 2 is regarded for by taking the previously determined values in the LED-delay experiment and in the image processing log data into account.

Results

The experiments on the feedback system are evaluated, starting from an elementary consideration of each element in the loop and its individual contribution to the feedback lag. Subsequently, the integral feedback lag is measured and its impact on the tracking experiment based on physiologic motion is shown.

Elementary timings of the feedback system

Qty. Table 1 shows the resulting lag contributions for each component individually. Very prominent is the relatively long line feed time of the imaging system ∆T 1 . Additionally, its distribution deviates significantly. The various tracking algorithms unsurprisingly yield different lags in the ∆T 2 -estimation. The correlation, being a global filter, is the slowest. The subscription delay ∆T 4 is extracted from a single parameter fit, thus no parameter statistics are available.

The open-loop delays ∆T 3 field are depicted in figure 4a in a map which represents the parameter field. Although, a Gaussian distribution cannot be presumed, overlapping median and mean ∆T 4 = ∆ T 4 = 15.9ms imply symmetry of the distribution with a estimated standard deviation of σ = 6.6ms. The minimal coefficient of determination of the obtained fits of the entire parameter field was R 2 min = 99.63%. The coefficients of determination indicate proper data fits and consequently significant lag extraction.

Closed-loop feedback lag determination

Figure 5 shows a section of the fitted data from object and aperture position. Both fits had R 2 > 99.94%. The extracted lag was determined to be ∆T ′ loop,corr = 109.3ms, based on the relative phase between the fitted curves and after correction according to equation 3. Irradiation under physiologic motion (blue traces in figure 2) yielded the normalized dose distribution in figure 6. As expected, the untracked distribution in Figure 6b shows significant smearing around the edges perpendicular to the leaf travel direction. Also lateral deviations are visible. Choosing an aperture relatively larger than the displacement of the tracked object (w aperture /2 > max{|∆x|}) yields the anticipated visible dose plateau phase in the center of the distribution. Figure 7 integrates the relative dose maps of figure 6 into histograms; 7a shows the distribution of dose difference between static object (reference) and untracked object. The multimodal distribution shows two lateral side maxima caused by the comparatively long excursion plateaus of the motion trace. In figure 7b, the distribution morphs into a compact, mean-free Gaussian-like distribution with σ (∆D) = 1.38% using tracking. In the γ-map of figure 8, The γ-pass rate, extracted from the thresholded (γ < 1) maps of figure 8 changed from 49% for the untracked case, to 92% for the tracked case. The concentration of γ-failed points is in regions of large spatial dose gradients, i.e. around the location of the leaf tips. The pass area in both γ evaluations at y = 0.02m shows intersecting dose profiles, the reference dose and the dose accumulating during tracking being at the same level. Due to the ideally strong gradients with the reference aperture a pass region will be observed with a width solely dependent on the geometric component of the γ-statistics.

Discussion

For the dynamic experiments, the Agility MLC was driven in a continuous manner, avoiding retarding inertial effects. This resulted in a low MLC delay of 15.9ms shown by indirect measurement. Although exhibiting significant variance, the phase featured linear phase behaviour in the sampled parameter field. Considering the real-time frame of 40ms, this small delay is likely to be caused by feed-forward terms in the MLC control system, imposing a ballistic extrapolation between adjascent control points of two consecutive cycles. These findings complement step-responses for the Varian Millenium MLC by (Poulsen et al. 2010) in which time delay was determined to be constant for a limited MLC adjustment step.

Aiming for potential MLC tracking in an MR-Linac, 2D imaging rates are expected to be in a range of 2-4 frames per second; accordingly an upsampled, extrapolated object position variable in order to update the MLC in every control system cycle could yield a dosimetric advantage of low MLC latency compared step-wise MLC response. Furthermore, using 1D navigator echoes, Stam et al. (2012) could get anatomy information every 15ms, which could be used as a fast surrogate tracking variable in combination with an organ correlation model.

Compared to the low apparent MLC latency, current clinical position extraction modalities such as approximately 90ms for Calypso (Smith et al. 2009), 100ms for MRI (Ries et al. 2010) and > 160ms for kV/MV imagers (Poulsen et al. 2010) contribute to the overall latency to a much higher extent. Notably, for the non-clinical tracking implementation with a customary camera, the most significant contributor to the overall delay of about 110ms was the imaging API itself. The API delivered images with an average of 78.7ms, which can be explained by the buffer strategy of the frame capturing cascade of the OpenCV API (Schellen et al. 2013). It is however expected that this lag component can be be overcome in a clinical modality using real-time communications as demonstrated by Ries et al. (2010) and Crijns, Raaymakers & Lagendijk (2012).

Knowing the MLC-inherent logging delay ∆T 4 in addition to the MLC delay, dose can be accumulated retrospectively by knowing both temporal and geometric interplay of MLC and anatomy. The variance of ∆T 4 , although being extracted from only a single indirect measurement, is expected to be small, since the position subscription runs with a time frame of the Agility MLC control system. Dosimetric deviations are distributed in a non-standard distribution for both tracking and static MLC applications. Qualitatively, both the relative differences and the γ-score show the performance improvement of tracking (γ = 92%) versus static (γ = 49%) MLC irradiation. Since the previously published, more relaxed γ statistics with radii of as shown by Smith et al. (2009) 3mm/3% are intended to rate the outcome of the entire treatment chain, the stringent γ (1mm/1%)-statistics may be more suitable as a representative integrated quantifier for rating the tracking performance of the tracking-system in an isolated manner. Observed dosimetric deviations laterally, perpendicular to the leaf travel direction, is due to to radial anisotropy of the beam itself and the dosimetry setup, i.e. differing path lengths through the dosimetry set-up. The latter is adding to the fact that retrospective dose accumulation cannot solely be computed with a simple convolution kernel constructed by system delay, but has to account for the beam characteristics, e.g. by incorporating the aperture and anatomy position into the treatment planning system (TPS) for retrospective dose calculation published by Poulsen et al. (2012) and clinically used by Keall et al. (2014). This is especially true in flattening filter-free applications. However, while knowledge about machine geometry is comparatively easy to obtain knowing the timing parameterization, more details about tissue position over time has to be integrated into the process either by fast 3D imaging techniques or precise motion models using surrogates. Recently, Bjerre et al. (2013) and Brix et al. (2014) could successfully employ MRI to rapidly infer 3D anatomy motion.

For applying tracking in the second dimension perpendicular to the leaf-motion, on-line sequencing needs to be elaborated. Also, the dependency of the MLC dynamics on gantry and collimator angle needs further investigation. As shown within this paper, position prediction filters will be needed to be further investigated in order to keep the MLC in the presented continous motion and to compensate for the potentially variable system delay.

Conclusion

The tested Agility MLC system showed notably low latency times when driven in a continuous manner; extracted timing characteristics indicates good MLC tracking performance with regards to commonly employed imaging and surrogate positioning techniques. Counteracting bottlenecks in tracking cascades will thus have to emphasize on fast and responsive imaging techniques. The shown phase linearity enables generic parameterization of predictive algorithms. An adaptive framework for dynamic multileaf collimator (dMLC) tracking employing the Agility MLC could succesfully be employed and first dosimetric experiments based on in-vivo liver motion showed expectably improved dosimetry performance both in relative dosimetry and in a γ-sense.

Figure 1 .

 1 Figure 1. Schematic (a) and picture (b) of system components of the feedback loop, information flow and motions: The controller is commanding the MLC system and retrieving leaf positions from it. The position extraction locates the object position by processing a video stream from a camera. 1D motion is synthesized and sent to the motion stage which translates the dosimetry stage and the attached target. See text for details.

  Figure 2. Breathing motion data of two volunteers and the (15mm, 6s) sinusoid in temporal (a) and frequency (b) domain. Notably, most signal energy is under 0.3Hz. The spectrum is scaled to fit the physiologic signals.

Figure 3 .

 3 Figure 3. System lag components and experimental pathways in the feedback and open-loop configuration. The upper solid path emanating from the detected object motion represents closed-loop (or tracking) experiments. The other, lower solid paths show open-loop experiments with directly controlled MLC. For determination of the logging delay, the dashed pathway was employed. Diamonds indicate position data sampled. δ M LC,logged and δ obj.pos denote logging timing for MLC and object, respectively.

Figure 4 .

 4 Figure 4. Open-loop MLC latencies for the sampled parameter field in (a). The distribution of the delay times is shown in (b).

Figure 5 .

 5 Figure5. Closed-loop logging data and relative lag determination using fitted objection (solid-blue with • data points) and aperture motion (red with ×-illustrated points. Solely every third data point is shown for clarity.

Figure 6 .

 6 Figure 6. Normalized dose distributions under physiologic conditions: in a static (a), untracked dynamic (b) and tracked dynamic case (c). Leaf travel direction is the vertical direction. Axis units in m.

Figure 7 .

 7 Figure 7. Relative difference histograms of untracked (a) and tracked (b) relative dose distribution.

Figure 8 .

 8 Figure 8. Thresholded γ (1%1mm) maps of films without (a) and with MLC tracking (b). The red area designates the region of D rel > 10% where the γ-test has passed, whereas white indicates points of failed test. Axis units in m.

Table 1 .

 1 Results for average values of dead times ∆T 1,2,3,4 with their estimated standard deviations σ.

		µ	σ	remark
	-	ms ms	-
	∆T 1 78.7 19.6	multimodal
	∆T 2 20.6 1.7	correlation (100 × 38px)
		5	0.2	MOSSE (32x32px)
		11.9 0.4	MOSSE (64x64px)
	∆T 3 15.9 6.6 sampled parameter field (figure 4)
	∆T 4 52	-	one experiment (param.-fit)

Acknowledgments

This work was supported by research funds from Elekta Ltd., UK. We appreciate their efforts in providing experimental research versions of machine firmware which enabled this work.


Juan Pablo Peralta

Maria Alejandra Aguirre

Jean-Christophe Géminard

Luis

Alejandra Aguirre

Luis A Pugnaloni
email: luis.pugnaloni@frlp.utn.edu.ar

Apparent mass during silo discharge: Nonlinear effects related to filling protocols

Keywords: latter puzzling experimental observation Granular flow, Granular discharge, Stress in silos

Introduction

Silos are very common structures that are broadly used to store powders or grains, i.e. granular materials. Because of the obvious importance of silos for the industry (e.g. pharmaceutical, mining, agriculture and others), questions related to the flow out of the structure through an aperture opened in the base and to the stress applied by the granular material on the walls of the structure have attracted the interest of engineers and physicists for decades.

On the one hand, when the granular material is discharged through an aperture in the base, one observes that the flow rate is constant throughout the discharge. Engineers and physicists have intensely studied the problem for the last 50 years [1][2][3][4][5][START_REF] Duran | Sands, Powders and Grains[END_REF][START_REF] Ristow | Pattern Formation in Granular Materials[END_REF][START_REF] Savage | [END_REF][START_REF] Hagen | Bericht uber die zur Bekanntmachung geeigneten Verhandlungen der Koniglich Preussischen Akademie der Wissenschaften zu Berlin[END_REF]. Heuristically, the flow rate, W , can be shown to scale as A 5/2 , where A is the size of the aperture [10][11][12][13]. This scaling, known since 1852 [START_REF] Hagen | Bericht uber die zur Bekanntmachung geeigneten Verhandlungen der Koniglich Preussischen Akademie der Wissenschaften zu Berlin[END_REF], is generally referred to as the Beverloo rule [1,[START_REF] Brown | Principles of Powder Mechanics[END_REF].

On the other hand, when granular material is stored in a vertical cylinder, the most remarkable feature is that the vertical stress applied by the material onto the base saturates to a constant value when the height of granular material becomes large compared to the diameter of the container. This is clearly in contrast with the pressure at the base of a viscous liquid column which increases linearly with the height of the liquid column. The effect, known as Janssen effect, is due to the fact that part of the weight of the grains is supported by the lateral walls of the container thanks to friction [START_REF] Hagen | Bericht uber die zur Bekanntmachung geeigneten Verhandlungen der Koniglich Preussischen Akademie der Wissenschaften zu Berlin[END_REF]15]. The Janssen effect has first been proposed to be at play in static conditions (in absence of flow) and reported by several authors in dedicated studies [16][17][18][19][20].

More recently, a screening of the weight of the grains at the bottom of a container has been observed in dynamic situations [21,22], suggesting that a dynamic Janssen-like effect can also be at play during a silo discharge.

One could thus infer that the flow-rate out of a silo is constant during the discharge, independent of the height of material above the outlet, because the pressure in the outlet region is constant [23]. However, measurements of the stress at the base of a discharging silo [21,[24][25][26][27][28] and numerical studies [29] showed that the flow-rate remains constant throughout the discharge in spite of significant variation of the pressure in the outlet region. Thus, unlike viscous fluid flows -which are controlled by the pressure-granular flows through an aperture are not sensitive to the local pressure.

Constitutive equations describing the flows of granular materials in all their generality are still lacking; as illustrated by the puzzling features exhibited by the old and, apparently, simple problem of the silo discharge. In order to provide additional clues to understand the problem, we consider the effects of the filling protocol which has been shown to affect the redirection of forces towards the walls [26]. In particular, we consider the evolution of the force exerted by the grains on the bottom surface of a discharging silo for systems prepared using a combination of two filling protocols. We will show that some unexpected nonlinear effects are observed when two filling procedures are used. However, the use of a modified Janssen model allows us to quantitatively explain the results.

Experimental setup and protocols

Experimental setup

The basic idea is to measure the overall vertical force exerted on the base of a vertical silo during the discharge of granular material prepared using various protocols.

The experimental setup is similar to the one decribed in [21]. It consists of a silo whose vertical outer wall is a Plexiglas cylinder of internal diameter D c = 5 cm, external diameter D e = 6 cm and total height H = 40 cm (Fig. 1).

At the bottom, an acrylic ring of internal diameter equal to D c insures the mechanical contact between the cylinder and a horizontal table, such that the vertical force applied to the outer wall (its weight plus the additional force exerted by the grains) is entirely supported by the table. Inside the latter outer ring, the base of the silo consists of an acrylic disk which leans on two force sensors (KD 40S-Testwell). The outer diameter of the disk is slightly smaller than the inner diameter of the outer ring, such that the two parts are not in mechanical contact, but the gap is small enough to avoid the grains to enter therein. In this configuration, the force sensors measure the total vertical force applied to the bottom, or equivalently the apparent mass, M app , of the material inside the silo. We shall report measurement of the apparent mass M app during the discharge of the silo. To do so, an orifice of diameter A = 5 mm has been drilled at the center of the base. In order to insure a sharp edge and avoid any effect of the thickness of the disk, the orifice is slightly conical, the small diameter in contact with the grains above.

During the filling procedure, the orifice is closed with a stopper which is held in place with a clamp.

The granular material consists of a total mass of 1 kg of glass beads (ρ glass = 2200 kg m -3 ) whose diameter, D, is broadly distributed between 300 µm and 600 µm.

The discharge is initiated by releasing the stopper, suddenly but without jerk thanks to the clamp. Underneath the aperture, the grains are collected in a container placed on an electronic scale. Thanks to the scale and force sensors, we measure simultaneously the discharged mass, M (t), and the apparent mass, M app (t), as function of time, t. In what follows, for each of the experimental protocols used to prepare the system, we report data obtained for 4 to 10 experimental runs. We observed that the results are reproducible for all the protocols. The error bars correspond to the standard error of the mean, which is calculated over the independent experimental runs.

Filling procedure

In the present study, we focus on the sensitivity of the apparent mass on the way the material is prepared previous to the discharge. We therefore pay special attention to the filling procedure. We will mainly report data obtained for silos prepared using 4 different ways, the total initial mass of material in the silo remaining the same, equal to 1 kg:

• Concentric (CC ): 1 kg of granular material is poured as a jet from the top along the vertical axis of the silo using a centered funnel;

• Distributed (DD): 1 kg of granular material is homogeneously deposited across the section of the silo by releasing the grains through three successive meshes [30];

• Mixed I (CD): Initially, half of the total mass of grains, i.e. 500 g, is introduced in the silo using the concentric filling procedure; the second half of the grains being introduced using the distributed filling.

• Mixed II (DC ): Initially, half of the total mass of grains, i.e. 500 g, is introduced in the silo using the distributed filling procedure; the second half of the grains being introduced using the concentric filling.

The filling procedures lead to different initial packing fractions, even if the system is filled at the same rate. For instance, the time required to fill the silo using the CC and DD procedures are the same for our choice of funnel and meshes. However, we measure that, in average, the initial height of the granular column is h CC ∼ = 33 cm for the CC filling whereas h DD ∼ = 30 cm for the DD filling. The observation is compatible with previous results reported in the literature [26] showing that DD filled systems are denser than CC filled systems.

Experimental results

In this section, we report on the evolution of M app during the discharge for systems filled with the different procedures described in Sec. 2.2.

We first comment that the flow-rate remains constant throughout the discharge and does not depend on the filling procedure. In addition, a convenient way to display the results consists in reporting M app as function of the mass of grains, M in , remaining inside the silo. The results are shown in Fig. 2. The initial mass of grains in the silo is M in = 1 kg in all cases. Therefore, the discharge proceeds from right to left in the figure. Notice also that the final mass of grains inside the silo does not vanish since a small quantity of grains remains in a pile around the orifice at the end of the discharge.

Single filling protocols

We first focus on the results obtained for systems prepared using a single filling procedure, i.e. concentric (CC ) or distributed (DD).

In the initial stages of the discharge (Fig. 2), M app is considerably lower for the CC than for the DD filled systems. At first sight, this observation contradicts results by Ovarlez et al. [20], where granular systems with a higher packing fraction present a larger redirection of the granular weight to the wall, i.e. smaller M app . As mentioned above, in our experimental conditions, the density achieved by the DD filling is larger than that achieved by the CC filling, which should lead to a smaller M app for the DD filling, in contradiction with our experimental observations. Note however that Ovarlez In contrast, we fill the silo using different protocols that lead to different densities of the granular packing, but also to different stress distributions. It is important to realize that the packing fraction is not an indicator of how the stress is distributed in the bulk of the sample and, therefore, can not be used as a proxy for how the weight of the grains is redirected toward the walls.

The different M app observed in the CC and DD filled systems can be qualitatively understood as follows. When the silo is filled using the concentric procedure, the grains initially form a pile at the center of the base.

Then, as grains are added to the summit, they produce avalanches that are stopped by the lateral walls. During the filling, this dynamics creates a rather large pressure on the silo walls. In contrast, when the silo is filled using the distributed procedure, the grains fall rather uniformly across the entire cross-section of the silo and the free surface remains flat throughout the filling. The grains do not induce avalances and one can speculate that they apply rather small normal forces on the silo walls. A large pressure on the wall leads to a large dynamic friction between the moving grains and the walls during the discharge. This qualitatively explains why the system prepared with the concentric filling has a larger portion of its weight sustained by wall friction, leading to a smaller M app , in comparison with the samples prepared using the distributed filling, even if the density of the material is smaller in the concentric filling case.

Mixed filling processes

Let us now discuss the results obtained for systems prepared using mixed filling procedures, i.e. concentric followed by distributed (CD) and distributed followed by concentric (DC ).

We focus first on the initial stage of the discharge in Fig. 2. We observe that the apparent mass M app is systematically larger when the lower half of the silo was filled using the distributed procedure, i.e. DD and DC, rather than using the concentric procedure, i.e. CC and CD. This observation can again be understood, as in Sec. 3.1, by considering that the redirection of the weight of the grains toward the lateral walls in th lower half is enhanced by the avalanching of the grains during the filling when the concentric (jet-like) protocol is used. We also observe that silos whose lower half is filled using the concentric procedure, i.e. CC and CD, present approximately the same M app at the very beginning of the discharge. This indicates that the redirection of the weight toward the silo walls is efficient and that the lower half of the silo redirects almost entirely the mass of the upper half toward the walls. In contrast, for silos whose lower half is filled using the distributed procedure, i.e. DD and DC, the initial apparent mass depends on the filling protocol used for the upper half. For the DC filling procedure, M app is smaller than for DD filling procedure. The observation is compatible with the fact that the apparent mass of the upper half supported by the lower half is smaller when the upper half is filled using the concentric procedure.

We now consider the evolution of the apparent mass throughout the discharge. We observe that the CC, DD, and DC fillings lead to monotonic dependencies of M app versus M in . However, the silo filled using the CD protocol has an non-monotonic behavior and M app exhibits a minimum which is considerably below the apparent mass measured for the CC filling. This experimental fact is rather unexpected since a naive interpolation argument would suggest that the CD filling procedure should lead to M app values in between the CC and the DD filling protocols. The result is particularly interesting when considering applications. Indeed, an unexpected small apparent mass is associated to an unexpected large pressure exerted on the lateral walls, which might be of practical importance for the design of silos.

Our experimental observations can be explained by using rather simple theoretical arguments. For instance, the very nonlinear response observed for the CD filling procedure can be explained by considering the internal rearrangements induced in the concentric lower half. However, before we go further in the analysis of the results, we propose additional experiments that provide clues to identify the relevant physical ingredients to understand the observations.

Additional experiments

In order to get further insight in the physical mechanisms that govern the apparent mass of the granular column during the discharge that is described in Sec. 3.2, we performed a series of complementary experiments that make possible to isolate effects like the interaction between the grains at the interface between the lower and upper half of the granular column, the friction at the wall, etc. More precisely, we analyzed the discharge of silos whose lower half was filled using either a concentric (C ) or a distributed (D) procedure, completed using one of the following:

• a block of metal, having a mass of 240 g on top. This mass is similar to the apparent mass exerted by 500 g of grains filled using the distributed procedure. We denote CW (DW ) for concentric (distributed ) filling of the lower half with the additional weight on top (Fig. 3b).

• reduced friction of the upper half with the lateral walls. To do so, after the lower half has been filled, a cylindrical acetate film is introduced in the silo (covering the lateral walls of the upper half) before the remaining 500 g of grains are poured in the system to fill the silo.

In this case, we denote as DC o a distributed filled lower half and a concentric filled upper half with reduced friction (Fig. 3c). The other configurations are denoted in a similar way.

• a cardboard disk separating the lower and upper halves filled with grains. We denote as D|C a distributed filled lower half and a concentric filled upper half separated by the cardboard disk (Fig. 3d). The other configurations are denoted in a similar way.

In the following sections 3.3.1 and 3.3.2, we shall report results obtained for, respectively, concentric and distributed filling of the lower half.

Lower-half filled using the concentric procedure

We show in Fig. 4 the evolution of M app as function of M in for systems prepared as decribed above with the lower half filled using the concentric procedure. The measurements are displayed assuming a total inner mass of 1 kg even when the metallic weight is used, leading all the curves to start from M in = 1000 g.

Let us fist focus on the initial stage of the discharge in Fig. 4a. We observe that the systems prepared with the solid block (CW ) and the ones with a distributed filling of the upper half [with (C|D) or without (CD) the cardboard separator] all lead to the same initial apparent mass M app once the discharge is initiated. This confirms our suggestion that, initially, the upperhalf filled with a distributed method acts as an effective overweight of about 240 g. After about 200 g of grains have been discharged, we observe in Fig. 4a that M app for CW and C|D preparation protocols remains higher than the apparent mass observed for the CD filling procedure. This demonstrates that preventing the direct interaction between the upper and lower half (by using the cardboard separator) has a noticeable effect.

When the wall friction is reduced in the upper part (Fig. 4b) by using the acetate film, the evolution of the apparent mass M app is similar to that observed with the solid block (CW in Fig. 4a). However, in these cases (i.e., CC o and CD o ), due to the reduction of the friction at the lateral walls, the apparent mass M app reaches a larger maximum value (of about 300 g) when the material of the lower half has been entirely discharged, i.e. for M in 500 g. It is worth mentioning here that, in Fig. 4a, M app for CW and C|D does not reach a maximum of 240 g due to the grains that remain in the silo after the discharge and that still partially redirect its weight toward the walls. Finally, let us comment that the concentric filling of the upperhalf (CC) leads to a larger apparent mass than the distributed filling (CD), which is quite surprising as the apparent mass of the upper half is expected to be smaller for CC. This indicates that the distribution of the stress at the interface between the lower and upper halves is such that the weight of the upper-half filled using the distributed procedure is somehow better redirected toward the lateral walls by the grains that are underneath.

Lower-half filled using the distributed procedure

In order to get further insight into the role played by the distribution of the stress at the interface between the lower and upper halves and by the friction at the lateral walls, we consider similar experiments with the lower half filled using the distributed procedure.

As can be observed in Fig. 5, even if the use of a distributed upper half (DD) leads to a larger apparent mass throughout the discharge, the use of the acetate film at the walls (DC o ) or of the cardboard separator (D|C) does not alter the apparent mass M app as much as previously observed for the concentric lower half (see Sec. 3.3.1). This conclusion also holds when the solid block is used (DW ). This may be due to the fact that the distributed filling of the lower part presents a rather flat free-surface when compared with the heap formed by the concentric filling.

While the reason for the lower M app in DC -compared to DD-is due to the upper concentric half being more favorably supported by the walls, in the case of DC o -where wall friction is nearly absent-the lower M app is caused by an effect similar to adding a solid overweight. This induces a compression of the upper layers of the lower half that leads to an increased pressure and hence higher frictional forces on the walls, reducing the apparent mass with respect to the DD case.

Analysis and model

As mentioned in Sec. 3, it is quite remarkable that we measure that, as observed in Fig. 2 once the discharge is initiated, a system fully filled using the concentric procedure (CC) leads to a larger apparent mass than a system consisting of a concentric filling of the lower half with the upper half filled using the distributed procedure (CD). This is however compatible with the additional observation that the redirection of the weight toward the walls, and thus the screening effect, is enhanced when the pressure on the lower half due to the upper half is larger (CD, CW , C|D and CD o in Fig. 4).

Qualitatively, one can infer that an increase of the pressure at the interface leads, in a thin layer of grains below the interface, to a significant increase of the radial stress applied to the lateral walls. Thus, an enhanced dynamic Janssen-like effect during the discharge is observed. As a consequence, a smaller M app is measured at the base. To model this effect, we propose a Janssen-like rationale where the different filling protocols lead to different redirection of the forces toward the lateral walls.

Janssen rationale

In a granular column, the pressure P (z) in the vertical direction, i.e. the normal component of the stress along the vertical axis, σ zz , is governed by the Janssen differential equation [15]:

∂P (z) ∂z -ρg + 2µk R P (z) = 0, ( 1
)
where z is the vertical position measured from the free surface downwards, g the acceleration of gravity, ρ the apparent density of the granular column, R the radius of the silo, µ the grain-wall friction coefficient and k the so called Janssen force-redirection factor. The latter constant factor k accounts for the proportion of the normal vertical stress transferred to the normal radial stress at the walls, i.e. σ rr = k σ zz . Usually, Eq. ( 1) is closed with the boundary condition P (0) = 0 (zero pressure at the free surface) and the well known Janssen law is obtained in the form

P (z) = ρgR 2µk 1 -e -2µkz R . (2)
From Eq. ( 2), the apparent mass M app = [πR 2 P (0)]/g as function of the

inner mass M in = ρπR 2 z is M app = M c 1 -e -M in Mc , (3)
where M c ≡ ρπR 3 /(2µk) is a characteristic mass.

As expected from Eq. ( 3), M app indeed saturates for large values of M in [see the fit of Eq. ( 3) to the CC data in Fig. 6]. The best fit to the experimental data to Eq. ( 3) for the concentric filling provides us with a measure of the characteristic mass, M c = (220 ± 5) g. From the knowledge of R = 0.025 m, g = 9.8 m/s 2 , ρ = 1320 kg/m 3 and M o = 1000 g, we estimate the product µ k.

Assuming further that µ = 0.51 , we get an estimate of the force-redirection factor k = 0.31 for the concentric filling (Table 1).

However, one immediately notice that Eq. ( 3) would fail in accounting for the dependence of M app on M in for the distributed filling procedure. Indeed, no saturation is observed (see DD in Fig. 6) and there is an almost linear increase of M app with M in .

Modified Janssen-like rationale for distributed filling

The only way to explain the linear increase of M app with M in for large inner mass is to consider that the redirection factor k almost vanishes for the distributed filling. However, some sort of screening is still at play as proven by the fact that the experimental slope of M app versus M in is about 1/3 (much less than unity). We thus propose that the redirection of the force toward the lateral walls is very inefficient far from the outlet, i.e. that the redirection factor k is almost zero along the column height for the distributed filling. However, the convergent flow of the granular material close to the outlet leads to a local redirection of the forces and, hence, the weight of the granular column is partly supported by the lateral walls in this region, only.

A convenient way to translate these ideas into a mathamtical model is to write that the redirection factor takes a finite value k 0 at the base and decreases over a characteristic vertical distance L above. We thus write

k(z) = k 0 exp(-z/L). ( 4
)
In this picture, the distance L is governed by the geometry and one expects it to be of the order of the container diameter 2R.

Following a Janssen-like rationale, we write, for the pressure along the vertical direction

∂P (z) ∂z -ρg + 2µk 0 R exp(-z/L) P (z) = 0, (5)
The solution of Eq. ( 5), with the additional boundary condition that P = 0 at the free surface, leads to

M app = M c Le -L E i L -E i L e - M in LMc , ( 6
)
where we introduced the reduced length L ≡ 2k 0 µL/R and the characteristic mass M c ≡ ρπR 3 /(2µk 0 ) defined using k 0 . The function E i (x) is the exponential integral given by

E i (x) ≡ - ∞ -x e -u
u du. Eq. ( 6) predicts, after a transient, the linear increase of M app as a function of M in with a slope smaller than unity. The best fit of the experimental data for the DD filling in Fig. 6 using Eq. ( 6) leads to estimates of both the redirection factor k 0 0.7 and the characteristic length L 4.4 cm (Table 1). Note that this latter value is, as expected if the model is correct, similar to the container diameter. Note also that the redirection factor k 0 corresponds to a reasonable angle of redirection of tan -1 (0.7) ≈ 35 deg (this is about the angle of repose of the pile that remains at rest around the outlet at the end of the discharge) Our Janssen-like rationale accounts thus correctly for the measurements of the apparent mass in the case of a simple distributed filling. In the next section, we show that the combination of a Janssen screening effect, with a constant value of k for the concentric filling (Sec. 4.1) and a "short-ranged" k(z) for the distributed filling (Sec. 4.2) accounts for the experimental measurements of the apparent mass when a mixed procedure is used.

Mixed filling procedure

We considered granular silos consisting of two halves, prepared in different ways. When solving Eq. (1) for a lower half filled using the concentric procedure or Eq. ( 5) for a lower half filled using the distributed procedure, one must thus consider as boundary condition, at the interface separating the two parts of the system, the pressure resulting from the upper part of the column of grains, instead of P = 0.

For a lower half filled using the concentric procedure, the solution of Eq. (1) with the condition that P (h lo ) = P 0 , where h lo denotes the height of the lower part of the column (which decreases during the discharge) and P 0 the pressure resulting from the grains of the upper half that lean on top, leads to

M app = M c 1 -e -2M in -Mo 2Mc + M up app e -2M in -Mo 2Mc . (7)
We have taken into account that the height, h lo of the lower part is a function of the inner mass according to ρπR 2 h lo = M in -M o /2 where M o is the initial total mass on grains poured in the silo. The relation holds true as long as M in > M o /2, i.e. while there remain grains of the lower half in the silo.

In Eq. ( 7), M up app is the apparent mass of the upper half estimated at the interface. For an upper half prepared with a distributed filling, M up app is given by Eq. ( 6) with M in = M o /2. This is because the total mass of grains in Table 1: Experimental redirection factor, k, (for concentrically filled part of the column) and characteristic redirection, k 0 , and characteristic length, L, (for the distributed filled part of the column) in the four different filling protocols.

(M in < M o /2) is then described by Eq. ( 6) when the grains of the upper part, filled using the distributed protocol, are flowing out of the silo. The best fit to the experimental data reported for the CD filling in Fig. 6 to Eqs. ( 6) and ( 7) leads to the redirection factor k 0.37 for the concentric half and to the redirection factor k 0 0.67 and the characteristic length L = 7.7 cm for the distributed half. These values, reported in Table 1, will be discussed below.

Finally, for a lower half filled using the distributed procedure, we can solve, in a similar way, Eq. ( 5) with the condition that P (h b ) = P 0 , where h b denotes the height of the lower part of the column, the pressure P 0 resulting from the apparent mass of the upper half, filled using the concentric protocol, which is given by Eq. ( 3) with M in = M o /2 during the first half of the discharge. The best fit of the experimental data reported for the DC filling in Fig. 6 to the model leads to the redirection factor k 0.21 for the concentric filling and to the redirection factor k 0 0.60 and the characteristic length L = 7.5 cm for the distributed filling (Table 1).

We have done an analogous analysis for the DC filling protocol, where the lower half is modeled using Eq. ( 5) with the boundary condition P 0 = M up app . The values obtained for k, K 0 and L for the best fit to DC data in Fig. 6 are listed in Table 1. We can see that, in spite of the experimental scatter, k, k 0 and L remain of the same order in all the experimental conditions used in our study.

We first comment on the values of the redirection factor, k, measured for the concentric filling. For the CC filling we measured k 0.31 whereas k is found slightly larger (k = 0.37) for a concentric bottom half underneath a distributed upper half and significantly smaller (k = 0.21) when the concentric half is on top of a distributed lower half. For the CD filling, the filling procedure is interrupted before the grains are poured as a rain on top. It might be that the mechanical perturbations associated to the manipulation have the effect of a slight tapping of the system, which might result in a more efficient redirection of the force toward the lateral walls (larger k), as it is observed when adding a solid overweight. For the DC filling, the discharge of the material underneath is likely to disorganize the texture of the concentric upper half, leading to a lower efficiency the redirection toward the walls (smaller k).

Let us now focus on the value of the redirection factor, k 0 , and of the characteristic length, L, measured for the distributed filling. We remind that we considered that the redirection of the force toward the lateral walls is induced by the geometry of the convergent flow at the lower region of the column of grains. For the DD and DC fillings, the pile of steady grains that surround the outlet might redirect the force toward the walls and we estimated from the measured values of k 0 that they correspond to a typical angle of about 35 deg, compatible with the typical angle of avalanche. We remark that for the CD filling, the concentric lower-half initially screens out the contribution of the distributed upper-half. Thus, the redirection factor k 0 for the distributed upper-half plays an important role only when the first half has been discharged and when the material of the upper-half enters in contact with the bottom. It is thus likely that, when the upper half filled using the distributed protocol reaches the bottom of the container, the same reasoning holds and the value of k 0 must be of the same order in this case.

We indeed measure that k 0 remains around 0.7 in all the configurations.

Concerning the characteristic length, the values of L for the DD and DC fillings are expected to be of the same order because the material of the bottom is confined in the same geometry. It is more surprising, at first sight, to get almost the same for the CD filling, for which the distributed upper-half sits on the lower half and not on the solid bottom. But, again, the redirection for the distributed upper-half is expected to play an important role only when the first half has been discharged and when the material of the upper-half enters in contact with the bottom. In this case, the geometrical conditions are the same. We notice, however, a factor of two in the characteristic length L, but, considering that the model is crude, the characteristic length remains of the order of the container diameter.

Discussion and conclusions

We assessed, during the discharge, the evolution of the apparent mass of the grains measured at the base of a vertical silo. We observed that M app versus M in depends drastically on the protocol used to fill the system. We considered silos filled in part using a concentric and in part using a distributed protocol.

In all cases, a large portion of the weight of the grains is supported by the lateral walls. We observe that the redirection of the weight toward the walls is more efficient for the concentric than for the distributed filling protocols.

For the concentric protocol, the redirection is so efficient that it leads to the saturation of the apparent mass when the height of the granular column is increased. In contrast, for the distributed protocol, the redirection is so inefficient that the apparent mass does not saturate when the height of the granular column is increased. From the theoretical analysis, we conclude that, for the distributed filling, about 1/3 of the weight is supported by the base of the silo, whereas the remaining part (about 2/3) of the weight is supported by friction over a small lower region of the wall whose height compares to the diameter of the silo. This implies that the grains exert a large normal stress onto the lateral walls only in this region for the distributed filling protocol.

We also found that filling the silo in two stages, half with concentric filling and half with distributed filling, leads to a nonlinear, non-monotonic evolution of M app as a function of M in . If the lower half corresponds to the concentric filling procedure, the system keeps memory of the filling history and shows a marked change in the pressure trend when the second half of material (filled with the distributed procedure) starts to evacuate. In addition, we observed that filling the container using two successive different protocols can lead to an apparent mass smaller than that achieved using any of the single filling protocols. This means that the lateral walls sustain larger stresses than those obtained with just one filling protocol.

We have shown that a Janssen-like model can be used to quantitatively explain the experimental results and to predict the apparent masses for the still little explored cases of distributed filling and mixed filling protocols.

These findings are of particular importance when designing operational protocols for filling and emptying silos, since the stress on walls can be largely affected and the system can keep memory of the history of the filling procedures even after significant manipulation. The larger the forces at the wall, the larger their wear and possible collapse [31]. Since it is generally simpler and more affordable to reinforce the base of a silo, a distributed filling protocol -which showed the higher M app -should be considered of preference over a concentric filling. However, one must consider to reinforce not only the bottom but also the lower part of the lateral walls in a region extending vertically over typically a silo diameter.

It is also important to mention that changes in the filling protocols must be considered with some care to avoid compromising a silo structure. If a silo usually filled concentrically is completed using a distributed procedurewithout emptying the remaining of previous concentrically filled materialthe walls may initially suffer a larger stress than keeping the usual concentric filling.

  et al. induce a relaxation of stresses to a steady state by slowly moving the bottom of the silo downwards by means of a vertical translational stage.

21 D

 21 the upper half remains constant and equal to M o /2 during the first half of the discharge (i.e. while M in > M o /2). The second half of the discharge CC DD CD DC C-filled part k 0.31 -0.37 0.

Figure 1 :

 1 Figure 1: Sketch of the experimental setup (see description in Sec. 2A).

Figure 2 :

 2 Figure 2: (Color online) Apparent mass M app vs. mass M in of grains inside the silo during the discharge for the different initial filling procedures described in Sec. 2.2: distributed (DD, red disks); concentric (CC, green squares); mixed I (CD, brown diamonds); mixed II (DC, blue triangles). The discharge proceeds from right to left. Error bars correspond to the standard error of the mean over 5 to 10 independent experimental runs.

Figure 3 :

 3 Figure 3: Additional experimental configurations -(a) Mixed II (DC ); (b) granular lower half plus a solid weight on top (DW ); (c) reduced friction of the upper half using an acetate film (DC o ); (d) lower and upper halves separated by a cardboard disk (D|C).

Figure 4 :

 4 Figure 4: (Color online) Apparent mass M app vs. inner mass M in for silos with the lower half filled using the concentric procedure -(a) concentric (CC, green squares); mixed I (CD, brown diamonds); concentric lower half plus overweight on top (CW, pink triangles); mixed I with a cardboard separator (C|D, orange disks). (b) concentric (CC, green squares); mixed I (CD, brown diamonds); concentric (CC o , pink triangles) and mixed I (CD o , orange disks) with reduced friction at the lateral walls in the upper half .

Figure 5 :

 5 Figure 5: (Color online) Apparent mass M app vs. inner mass M in for lower half filled using the distributed procedure -distributed (DD, red disks); mixed II (DC, blue triangles); mixed II with reduced friction at the lateral walls in the upper half (CD o , pink squares); mixed II with a cardboard separator (D|C, orange diamonds); distributed lower half plus overweight on top (DW, green crosses).

Figure 6 :

 6 Figure 6: (Color online) Apparent mass M app vs. inner mass M in . The symbols are as in Fig. 2, solid lines correspond to the fits of the modified Janssen-like model.

The value of µ is rather arbitrary. However, we use the same value in all cases (all the experiments where done with the same materials) simply to focus our discussion on k.

Acknowlegments -We thank M. E. Fernández for help with some of the experiments. This work has been supported by the ANPCyT (Argentina) through grant PICT-2012-2155, and the Franco-Argentinian International Associated Laboratory 'Physics and Mechanics of Fluids' CONICET-CNRS.

[1] W. A. Beverloo, H. A. Leninger and J. van de Valde, Chem. Eng. Sci.


Cyril Lorenzato

Chris Oerlemans

Alexandru Cernicanu

Mario Ries

Baudouin Denis De Senneville

Chrit Moonen

Clemens Bos
email: c.bos@umcutrecht.nl

Rapid dynamic R1/R2*/temperature assessment: a method with potential for monitoring drug delivery

Keywords: Rapid Dynamic Temperature/R 1 /R 2 Magnetic resonance imaging (MRI), apparent proton density (APD), drug delivery systems (DDS), thermosensitive liposome (TSL), relaxometry, R 1 , R 2 *, focused ultrasound, High Intensity focused ultrasound (HIFU

INTRODUCTION

Undesired side effects associated with chemotherapeutic drugs can be substantial.

Therefore, there is a need for drug delivery systems (DDS) that may favorably alter the pharmacokinetics of the drug [START_REF] Allen | Drug Delivery Systems: Entering the Mainstream[END_REF] and reduce the toxic side effects as observed with systemic chemotherapy. Taking advantage of the ability of DDS to carry both imaging and therapeutic agents, detectable drug carriers have been prepared that allow imaging the drug distribution non-invasively [START_REF] Xie | Nanoparticle-based theranostic agents[END_REF][START_REF] Lammers | Theranostic Nanomedicine[END_REF][START_REF] Kelkar | Theranostics: Combining Imaging and Therapy[END_REF]. In this context MRI is a potentially interesting modality: DDS that behave as temperature-sensitive switchable MR contrast agents have been prepared, and would allow monitoring their availability [START_REF] Terreno | Paramagnetic liposomes as innovative contrast agents for magnetic resonance (MR) molecular imaging applications[END_REF][START_REF] Lorenzato | MRI contrast variation of thermosensitive magnetoliposomes triggered by focused ultrasound: a tool for image-guided local drug delivery[END_REF] and their release process [START_REF] Fossheim | Thermosensitive paramagnetic liposomes for temperature control during MR imaging-guided hyperthermia: In vitro feasibility studies[END_REF][START_REF] Smet | Temperature-sensitive liposomes for doxorubicin delivery under MRI guidance[END_REF][START_REF] Zwart | On the feasibility of local drug delivery using thermo-sensitive liposomes and MR-guided focused ultrasound[END_REF] based on MR-relaxometry measurement of the transversal relaxation rates (R 2 , R 2 * ) and the longitudinal relaxation rates (R 1 ), respectively.

In addition, MRI can non-invasively provide temperature maps, e.g. using the Proton Resonance Frequency Shift method [START_REF] Poorter | Noninvasive MRI Thermometry with the Proton Resonance Frequency (PRF) Method: In Vivo Results in Human Muscle[END_REF], and can thus be used for continuous guidance of the hyperthermia procedure that triggers drug release [START_REF] Quesson | Magnetic resonance temperature imaging for guidance of thermotherapy[END_REF]. One of the techniques to achieve local hyperthermia is High Intensity Focused Ultrasound (HIFU), which has been successfully integrated with MRI, such that MR thermometry controls the HIFU heating (12). Therefore, it would be of great value for the development of these image guided drug delivery techniques to have access to quantitative relaxometry [START_REF] Fossheim | Thermosensitive paramagnetic liposomes for temperature control during MR imaging-guided hyperthermia: In vitro feasibility studies[END_REF][START_REF] Bos | Simultaneous monitoring of temperature and T1: Methods and preliminary results of application to drug delivery using thermosensitive liposomes[END_REF] and temperature at the same time.

In this work, a fast dynamic MRI method to simultaneously acquire PRF-based thermometry data, R 2 * and R 1 maps using a dynamic multi gradient echo sequence is presented. To assess the feasibility of this approach, the signal change in response to changes of R 1 and the influence of the variance of transverse relaxation rate R 2 * estimation on the quantification of the longitudinal relaxation rate R 1 were analyzed. Subsequently the method was validated in vitro in a porcine muscle sample, initially using a water bath for heating and then applying local heating using MR-guided HIFU. Finally, the method was used to monitor the release of a paramagnetic contrast agent from thermosensitive liposomes (TSL). The method is a step towards monitoring two important phases of imageguided drug delivery: availability of nanocarriers and release in the target tissue.

METHODS

Data Acquisition

Real time mapping of relaxation times R 1 and R 2 * was based on a multi-echo RFspoiled gradient echo sequence (ME-SPGR). For this sequence, the signal is given by, where S n (TE) is the signal of dynamic n for a given echo time TE, α is the flip angle, M 0 is the nuclear magnetization which can be thought of as an apparent proton density (APD), γ the gyromagnetic ratio, and TR the repetition time. For slice selective sequences, the sincterm reflects the effect of the macroscopic main field inhomogeneity, which is assumed for these anisotropic voxels to be dominated by the field variation ∆B 0 /dz in slice direction [START_REF] Dahnke | Limits of detection of SPIO at 3.0 T using T2* relaxometry[END_REF].

With a fixed TR and α, the ME-SPGR does not allow quantifying R 1 . Instead, a baseline R 1 map is measured first, which is then updated based on the signal changes for each ME-SPGR dynamic. In this study, the initial R 1 map was measured using a Look-Locker acquisition (15) (resolution = 1 × 1 × 3 mm 3 ; FOV = 128 mm; TR = 4s; α = 5°; Images were acquired for 30 time points after the inversion, starting at 30 ms with 60 ms intervals; TFE factor = 5).

Subsequently, a B 0 map was estimated from the phase images of a multislice ME-SPGR scan without slice gap. Phase evolution from echo to echo was limited to be smaller than  ans echoes with signal magnitude below than 10% of the first echo were excluded as proposed in [START_REF] Dahnke | Limits of detection of SPIO at 3.0 T using T2* relaxometry[END_REF]. The resonance frequency was then obtained by calculating the slope of a linear fit to the signal phase as a function of echo time. Next, A map of the local gradient field ∆B 0 /dz was derived from a linear fit to the frequency values in the slice direction.

Finally, the central slice of the ME-SPGR was dynamically acquired to map temperature, R

Data analysis and Processing

Dynamic R * 2 mapping R 2 * maps were calculated from the multi-echo data using Eq. [START_REF] Allen | Drug Delivery Systems: Entering the Mainstream[END_REF], assuming fixed α, M 0 and R 1 values: , Here, S 0 is the signal when TE approaches 0. Three ways of handling ∆B 0 /dz compensation for fitting R 2 * were evaluated: 1) ∆B 0 /dz is assumed to be zero and R 2 * can be extracted from a mono-exponential Levenberg-Marquardt (LM) fit; 2) ∆B 0 /dz is assumed to be constant; it is measured once using the multislice ME-SPGR phase maps and then the sinc term from the Eq. [START_REF] Xie | Nanoparticle-based theranostic agents[END_REF] is calculated and included in the R 2 * fitting process; 3) ∆B 0 /dz is updated for each dynamic using the method described by Schaeffter et al. [START_REF] Dahnke | Limits of detection of SPIO at 3.0 T using T2* relaxometry[END_REF]. In this last case, the measured ∆B 0 /dz map was used as an initial value for an iterative algorithm that updates the ∆B 0 /dz map from frame to frame. Briefly, the magnitude intensity for each echo is corrected (divided) by the corresponding absolute sinus cardinal term calculated from the first ∆B 0 /dz map. Then a first R 2 * estimate is obtained by fitting the ∆B 0 /dz corrected signal to the remaining mono-exponential decay term in equation [START_REF] Xie | Nanoparticle-based theranostic agents[END_REF]. Finally, based on this first R 2 * estimate, ∆B 0 /dz is adjusted and the process is iterated to minimize the squared error in [START_REF] Xie | Nanoparticle-based theranostic agents[END_REF].

Dynamic R1 mapping

The implemented method estimates a change in R 1 based on the signal magnitude change under SPGR conditions. To allow measuring absolute R 1 values, an initial R 1 map is needed, which was calculated from the Look-Locker data [START_REF] Look | Time Saving in Measurement of NMR and EPR Relaxation Times[END_REF]. R 1 change (∆R 1 ) over time was calculated using the signal change at the shortest TE. For that purpose the signal from each new dynamic is divided by the signal of the reference dynamic S 0 .

From Eq. [START_REF] Lammers | Theranostic Nanomedicine[END_REF] it is apparent that for quantitative R 1 measurements, influence of temperature change and R 2 * change should be considered. The APD is proportional to the equilibrium nuclear magnetization M 0 , with M 0 (T) = Nγ 2 h 2 I(I+1)H 0 / 3κT, where κ is Boltzmann's constant, N is the number of nuclear spins per unit volume, and T is the absolute temperature. Thus, M 0 is inversely proportional to the temperature, with a resulting signal change of approximately 0.29%/°C in the experimental range of 20 to 60 °C [START_REF] Chen | Investigation of proton density for measuring tissue temperature[END_REF].Hence, the APD influence on the signal ratio can be easily compensated if the relative temperature variation is known. The second factor in Eq. [START_REF] Lammers | Theranostic Nanomedicine[END_REF] is related to R 2 * and can be calculated from the dynamic R 2 * mapping. The ∆R 1 is then determined numerically, in an iterative fashion, from Eq. [START_REF] Lammers | Theranostic Nanomedicine[END_REF], which can be rewritten as follows after compensation for M 0 and R 2 * :

Using R 10 as obtained from the Look Locker method, this provides an R1 value at every time point.

Proton resonance frequency based thermometry

PRF shift was calculated using the unwrapped phase change, ΔΦ, at the last echo time of the dynamic ME-SPGR data. The temperature change ΔT n is related to ΔΦ by: [START_REF] Poorter | Noninvasive MRI Thermometry with the Proton Resonance Frequency (PRF) Method: In Vivo Results in Human Muscle[END_REF], where β the temperature dependence water chemical shift of -0.0101 ± 0.0004 ppm/°C [START_REF] Peters | Ex vivo tissue-type independence in protonresonance frequency shift MR thermometry[END_REF]. PRF-based thermometry was corrected for potential B 0 drift by using the phase of an agar gel (2% w/w R 2 * ≈ 20 s -1 ) that was at constant temperature as a reference.

ΔT n = ΔΦ / β•γ•TE•B 0

Efficiency assessment of the proposed method

Sensitivity of signal to R1 changes

The signal change in response to R 1 -changes was evaluated to assess the sensitivity of the method. To this end, assuming constant R 2 * , realistically detectable magnitude signal variations from -20% to 20%, were evaluated for different R 1 starting values (ranging from 1 to 5 s -1 ). Eq. [START_REF] Kelkar | Theranostics: Combining Imaging and Therapy[END_REF] was solved iteratively to determine the required R 1 change.

Influence of ∆R2* on R1

Eq. [START_REF] Kelkar | Theranostics: Combining Imaging and Therapy[END_REF] implies that there is a correlation between errors in the R 2 * -change estimate and the estimate of R 1 . For this we may quantify the influence of ΔR 2 * as R 1 can be solved from Eq. [START_REF] Kelkar | Theranostics: Combining Imaging and Therapy[END_REF] (see supplementary material):

from which the partial derivative can then be calculated.

Porcine muscle heating

A porcine muscle sample of 3 x 3 x 1 cm 3 was heated from 26 to 34ºC in a water bath. During heating, real time ME-SPGR images were acquired for 20 minutes. Looklocker based R 1 maps were calculated before heating and after the sample reached thermal equilibrium in order to validate the real time measurements. The variation of R 2 * and R 1 over time were determined in an ROI encompassing the sample. A Luxtron optical probe was placed on top of the sample to monitor temperature (Fig. 1b). The temperature dependence of the relaxation rates R 1 and R 2 * of the porcine muscle sample were determined using weighted linear regression of the relaxation rate to the temperature as measured by the PRFS method.

MR-guided HIFU heating

Relaxometry and thermometry in the range of mild hyperthermia in porcine muscle tissue were finally demonstrated using a clinical MR-HIFU therapy system (Sonalleve;

Philips Healthcare). HIFU heating was performed with an acoustic power of 20W for 30 s, starting after 30 s of scanning which would typicaly lead to a temperature rise of 20°C.

Beam-steering of the focused acoustic field was performed to create a 4-mm diameter heating cell [START_REF] Partanen | Reduction of peak acoustic pressure and shaping of heated region by use of multifoci sonications in MR-guided high-intensity focused ultrasound mediated mild hyperthermia[END_REF]. The mean and the standard deviation, from a ROI of 4x4 voxels centered on the heating pattern, were used for monitoring R 1 , R 2 * and temperature over time. LLbased R 1 mapping was performed to validate R 1 values at the end of the procedure. In addition, precision of R 2 * and R 1 for muscle was estimated before heating in the same region of the sample from the standard deviation over time for the cases when ∆B 0 /dz was assumed to be 0, assumed to be constant , or continuously updated, respectively. A fiber optic probe provided the baseline temperature measurement.

Liposome preparation

Gadolinium-containing thermosensitive liposomes (Gd-TSL) were prepared using the conventional thin-film hydration technique as described previously [START_REF] Koning | Selective transfer of a lipophilic prodrug of 5fluorodeoxyuridine from immunoliposomes to colon cancer cells[END_REF]. The resulting liposome dispersions were sized with sequential extrusion using a Lipex Extruder (Northern Lipids Inc., Vancouver, Canada) and polycarbonate membrane filters (Poretics Corporation, Livermore, CA) with pore diameters of 600, 400 and 200nm. Extra-liposomal ProHance was substituted by HBS using PD-10 columns (GE Healthcare, Little Chalfont, UK) for 3 consecutive times.

Liposome characterization

The average hydrodynamic size and polydispersity index (PDI) of the Gd-TSL was determined with dynamic light scattering (DLS) using a Malvern ALV CGS-3 system (Malvern Instruments Ltd., Worcestershire, United Kingdom).

The chelated Gd concentration was estimated from the measured relaxivity of the free Gd.

The r 1 relaxivity of the chelated gadolinium was calculated from the linear regression of the R 1 measured for different concentrations (ranging from 0.1 to 1 mM). Then the stock Gd-TSL solution was diluted 10 times in 3 samples of 2 ml. Each sample was heated for 10 minutes at 45ºC in presence of Triton X-100. R 1 was measured for each sample and both the concentration and the encapsulation rate were calculated.

Gd-HPDO3A Release from Thermo-Sensitive Liposomes

The feasibility to monitor the release of a paramagnetic contrast agent from thermosensitive liposomes (Gd-TSL), using this method was investigated. For that purpose two tubes were placed in a MR compatible water bath and heated from 38.6 to 46ºC during 20 minutes of scanning. One tube contained 1 mL of the Gd-TSL preparation diluted 10 times, corresponding to a total [Gd] = 0.4 mM, preheated to 60ºC for 15 minutes and served as a control. The other contained 1 mL of fresh Gd-TSL. The variation of the mean and standard deviation of the relaxation times R 1 and R 2 * and the PRF-based temperature was analyzed using two ROIs of 36 voxels (see Fig. 1c). A fiber optic probe in between the tubes provided independent temperature measurement.

Porcine muscle heating

The water bath study was performed to validate the real-time R 1 measurements (Fig. 3). The mean values of the LL-based R 1 measurements were 1.26 ± 0.04 s -1 and 1.17 ± 0.04 s -1 at 26℃ and at the thermal equilibrium of 34℃, respectively. With an SNR of 38 on the first echo of the ME-SPGR data, the near real time R 1 measurement after heating was 1.18 ± 0.07 s -1 and showed a good correspondence to the LL-based R 1 measurement. R 1 was found to decrease with a temperature dependence of -0.0161 ± 0.0005 s -1 .℃ -1 which corresponds to a change in T 1 of 1.31 ± 0.05 %.℃ -1 .

The value of R 2 * , when ∆B 0/dz was assumed to be constant, was on average 0.92 ± 0.16 s -1 lower than if ∆B 0 /dz was assumed to be zero, confirming the contribution of the field gradient to R 2 * .(Fig. 3b). However, without ∆B 0 /dz compensation, the measured value was more precise, with a standard deviation of 0.28 ± 0.17 s -1 vs 0.66 ± 0.21 s -1 with compensation. The computation time of the R 1 map, the R 2 * map and the temperature map was under 0.8s per frame allowing for near-real-time monitoring.

MR-guided HIFU heating

The areas with R 1 and R 2 * changes, and the heated area based on the PRF values showed a good correspondence in position and shape (Fig. 4 and Table 1). When taking into account APD changes with temperature, the bias in R 1 was essentially removed. The precision of R 1 was essentially independent of the way R 2 * was calculated, which may be related to the high SNR of 85 in this experiment, as well as the moderate R 2 of muscle.

Nevertheless, iteratively updating ∆B 0 /dz led to a higher variance in the R 2 * estimates than assuming constant or zero ∆B 0 /dz. PRF-based temperature precision was 0.37 ± 0.2°C. The HIFU heating of 30 seconds caused a temperature increase of 20ºC. Both, R 2 * and the susceptibility corrected R 2 * increased during the heating period and decreased during the cooling process. The final R 1 measurement from the Look-Locker, corresponded well to the last incrementally calculated R 1 determined from the dynamic series. Similar temperature dependence as measured in the water bath experiment were found for muscle: -0.0127 ± 0.0003s -1 .°C -1 and 0.35 ± 0.01s -1 .°C -1 for R 1 and R 2 * , respectively.

Gd-HPDO3A loaded temperature-sensitive liposomes

Gd-TSL encapsulating 60 mM [Gd(HPDO3A)(H2O)] were prepared using the lipid film hydration technique followed by sequential extrusion. The hydrodynamic diameter of the Gd-TSLs was 194 ± 18 nm with a corresponding PDI of 0.04. The melting phase transition temperature was 42.2 °C with a DSC peak maximum of 44.3 °C. The R 1 relaxation rate of the released gadolinium from the 10 times diluted stock Gd-TSL was on average equal to 3.1 ± 0.1 s -1 corresponding to a concentration of 4.0 ± 0.2 mM of the stock solution and an encapsulation rate of 6.6%.

Gd-HPDO3A Release from Thermo-Sensitive Liposomes

Finally, it was demonstrated that the proposed method can be used to monitor the release of a Gd-chelate from TSLs during a water bath heating experiment (Fig. 5). The need for the APD correction was underlined by the comparison of the R 1 calculated from the ME-SPGR data during heating, 2.01 ± 0.10 s -1 with vs 1.93 ± 0.11 s -1 without APD correction, and the Look Locker data acquired at the end, 1.99 ± 0.06 s-1. PRF-based thermometry from the TSL tube and the fiber optic thermometer measurements showed an r 2 coefficient of 0.98, validating the temperature measurements using this method.

For the freshly prepared Gd-TSL, a positive R 1 variation was observed from 38.6ºC to 42ºC (t = 11 min). During the phase transition, from 42.2ºC to 44ºC, a sigmoidal-like increase of R 1 was measured. Above the phase transition, R 1 values were similar to these of the preheated Gd-TSL control sample until the end of the experiment. The control sample presented a small negative temperature dependence, during the whole heating process.

Corresponding R 2 * measurements (Fig. 5e,f) a detecteble but small change was found upon release, which is probably related to the limited R 2 * effect of chelated Gd.

for correction of the bias that was present without taking APD changes into account in the Gd-TSL release experiments. However, this method requires a reference R 1 map to allow calculating absolute R 1 values from the magnitude signal changes, as well as a phase reference for the temperature updates. This introduces motion sensitivity, and may thus limit its applications as this method will suffer from the defects as the PRF method.

Moreover, Hijnen et al. have shown that change of Gd concentration biases the PRF thermometry mostly due to changes of bulk magnetic susceptibilities [START_REF] Hijnen | The magnetic susceptibility effect of gadolinium-based contrast agents on PRFS-based MR thermometry during thermal interventions[END_REF] similarly to the motion sensitivity of the PRFS method that is based on phase subtraction. By monitoring R 2 * over time, we may detect when the Gd concentration is stable, and then initiate the heating under control of thermometry.

With respect to the microscopic local field inhomogeneities, the R 2 * changes, albeit weak, (Fig. 5e,f) and did not exclude accurate PRF-thermometry and R 1 measurement.

Measured R 1 temperature dependences of ex vivo porcine muscle were found to be independent of the heating device (HIFU vs water bath heating). Comparable T 1 temperature dependence , 1.4 ± 0.2%.℃ -1 , were found in literature [START_REF] Cline | MR temperature mapping of focused ultrasound surgery[END_REF] even if such values are highly dependent on measurement technique and tissue composition [START_REF] Bottomley | A review of normal tissue hydrogen NMR relaxation times and relaxation mechanisms from 1-100 MHz: dependence on tissue type, NMR frequency, temperature, species, excision, and age[END_REF]. The computation time was substantially smaller than the heating and release processes, as well as the dynamic scan time of the multi-echo SPGR sequence and thus suited for near-realtime monitoring.

Properly quantifying the released fraction of a contrast agent from thermosensitive nanocarriers looking at R 1 variations is difficult. Indeed, these variations also reflect the temperature and the concentration changes which are tissue dependent. However, looking to R 2 * variation should allow us to look not only to the ΔR 1 but also to the ratio ΔR 2 * /ΔR 1 as proposed by Terreno et al. to measure e.g. pH using Gadolinium complexes [START_REF] Aime | A R2/R1 ratiometric procedure for a concentration-independent, pH-responsive, Gd(III)-based MRI agent[END_REF].

CONCLUSION

We presented and validated a method for simultaneous PRF-based thermometry, and R 1 and R 2 * mapping with a spatial and temporal resolution high enough to characterize the release of contrast agents from TSL and to monitor HIFU heating, with potential for application in MRI-guided drug delivery studies. Computation time was < 0.05 ms.voxel -1 witch allow for a complet mapping of a FOV equal to 128 x 128 in less than a second.

Figure captions

  The TSL formulations consisted of DPPC, DSPC, cholesterol and PEG2000-DSPE in a molar ratio of 67:15:13:5. The lipid concentration was 40 mM. The lipid mixture was dissolved in ethanol (10 mL) and evaporated to dryness by rotary evaporation under vacuum (Rotavapor R-210, BUCHI Laboratory Equipment, Zurich, Switzerland) and the resulting lipid film was further dried under a stream of N 2 . A solution containing 60 mM of Prohance (gadoteriol, Bracco Diagnostics Inc. NJ, USA) and 10 mM of HEPES Buffered Saline (HBS), containing 135 mM NaCl, pH adjusted to 7.4 was used to hydrate the lipid film.

Figure 1 :

 1 Figure 1: Schematic of MR scanning protocol (a) with the corresponding temperature vs. time profile for the water bath and the HIFU heating studies. Magnitude images of the experimental setup for water bath heating of the porcine muscle sample (b) and the TSL (c).

Figure 2 :Figure 3 :

 23 Figure 2: Simulated magnitude signal variation ΔSn/S0 and the corresponding increase/decrease of R1 for different R 1 starting values (color-coded scale) (a). Influence of uncertainties in the estimation of R 2 * change on the R 1 estimation (b), for an example starting value of R 1 = 1.25 s-1, expressed as the partial derivative ∂ R 1 /∂ ΔR 2 * for a TR = 50 ms, TE = 4.6 ms

Figure 4 : 30 Figure 5 :

 4305 Figure 4: Dynamic magnitude image measured at 1.5 minutes with ΔR 1 overlay (a) and ΔR 2 * overlay (b). The inset in (a) shows PRF thermometry data. Evolution over time of R 1 (c) and R 2 * (d), respectively, measured at the focal point (4x4 mm 2 ) during the HIFU heating applied from 0.5 to 1 min.

TableTable 1 :

 1 Comparison of ending R 1 values and the Look Locker measurements from the HIFU experiment

	Handling of ∆B 0/dz	Ending R 1 ] M 0 = constant 1 [s -	Bias a	Ending R M 0 (Temperature) 1 [s -1 ]	Bias b	R 1 Precision [s -1 ]	c	* R 2 Precision [s -1 ]	c
	ΔB	0 /dz = 0	1.538 ± 0.081 -0.0072	1.548 ± 0.08	0.0025 0.039 ± 0.007	0.725 ± 0.170
	ΔB	0 /dz constant	1.536 ± 0.082 -0.0079	1.547 ± 0.08	-0.0003 0.039 ± 0.007	0.724 ± 0.168
	ΔB	0 /dz updated	1.538 ± 0.082 -0.0077	1.55 ± 0.08	0.0003 0.042 ± 0.008	0.863 ± 0.198

a Bias = ending R 1 with M 0 constant -Look Locker R 1 measurement. b Bias = final R 1 with M 0 (T) -Look Locker final R 1 measurement. c Temporal standard deviation before heating.

ACKNOWLEDGEMENTS

This work was supported by the Sound Pharma project (funded by the European Research Council). We would like to extend our profound and warmest thanks to the regretted Dr.

Pierre Smirnov without whom this work would not have been possible.

RESULTS

Sensitivity of signal to R 1 changes and influence of R 2 * on R 1

The imposed variation of S n /S 0 and the corresponding resulting R 1n values are shown in Fig. 2a. For the scan parameters used in the experiments and for typical R 1 (0.8 to 2.0 s -1 ), the calculated ΔR 1 is almost directly proportional to the signal variation.

Consequently, assuming known R 2 *

, with e.g. R 1 = 1 s -1 (T 1 = 1000 ms) and a signal change of ± 3% (SNR ≈ 33), the detection limit for R 1 changes are around ± 0.04 s -1 (± 42 ms).

The partial derivative ∂R 1 /∂∆R 2 * was calculated (Eq. [START_REF] Terreno | Paramagnetic liposomes as innovative contrast agents for magnetic resonance (MR) molecular imaging applications[END_REF]) for in-vivo typical R 2 * values (from 10 to 60 s-1 = 100ms) to quantify the effect of errors in R 2 * estimation on the R 1 estimation (Fig. 2b). On one side, high R 2 * values > 50 s -1 lead to high ∂R 1 /∂R 2

DISCUSSION

In this study, it has been shown that R 1 , R 2 * and temperature can be monitored based on a multi-echo SPGR sequence.

The ability of a T 2 * shortening agent to create contrast even when encapsulated into nanoparticles [START_REF] Lorenzato | MRI contrast variation of thermosensitive magnetoliposomes triggered by focused ultrasound: a tool for image-guided local drug delivery[END_REF][START_REF] Bulte | Magnetoliposomes as contrast agents[END_REF][START_REF] Girard | Toward absolute quantification of iron oxide nanoparticles as well as cell internalized fraction using multiparametric MRI[END_REF][START_REF] Billotey | Cell internalization of anionic maghemite nanoparticles: Quantitative effect on magnetic resonance imaging[END_REF], makes R 2 * an interesting mechanism for monitoring the availability of drug-nanocarriers in the tissue. In addition, earlier studies using thermosensitive liposomes have shown the potential of MRI to monitor the release from nanocarriers using R 1 mapping [START_REF] Smet | Temperature-sensitive liposomes for doxorubicin delivery under MRI guidance[END_REF][START_REF] Viglianti | Chemodosimetry of in vivo tumor liposomal drug concentration using MRI[END_REF][START_REF] Ponce | Magnetic resonance imaging of temperature-sensitive liposome release: drug dose painting and antitumor effects[END_REF]. Dahnke et al. showed that for slice selective multi-echo sequences R 2 * has to be compensated for B 0 gradients in slice direction [START_REF] Dahnke | Limits of detection of SPIO at 3.0 T using T2* relaxometry[END_REF]. In our case, ∆B 0 /dz was weak and stable for temperatures in the mild hyperthermia range up to 43°C. However, as demonstrated in this paper, R 2 * changes need to be estimated accurately and precisely to allow for R 1 monitoring based on signal changes. In this study, three ways of handling B 0 in R 2 * estimation have been compared, and it was confirmed that including ∆B 0 /dz in the R 2 * estimation removes biases. Nevertheless, it may be advantageous to measure ∆B 0 /dz once and then keep it fixed, rather than to update it over time. Firstly, this allows using a single slice sequence for the dynamic R 2 * measurement. Secondly, in situations with low SNR, attempts to fit ∆B 0 /dz from the signal magnitude evolution with TE in a single slice for each dynamic will increase the variance on the R 2 * estimate, which may cross over to the R 1 estimates as well.

The PRF-based temperature measurements in turn allowed correcting APD changes in the R 1 estimate during heating. Obeying the Boltzmann distribution, the dependence of APD on temperature is well-known, and within the examined temperature range effectively allowed