Ph.D Jonathan Vappou
email: jvappou@unistra.fr

Pierre Bour

Fabrice Marquet

Valery Ozenne

Bruno Quesson

Quantitative MR-ARFI for HIFU therapy monitoring MR-ARFI-based method for quantitative measurement of tissue elasticity: Application to monitoring of HIFU therapy

Monitoring thermal therapies thanks to medical imaging is essential in order to ensure that they are safe, efficient and reliable. In this paper, we propose a new approach, halfway between MR Acoustic Radiation Force Imaging (MR-ARFI) and MR Elastography (MRE), allowing for quantitative measurement of tissue's elastic modulus in a highly localized manner. It relies on the simulation of the MR-ARFI profile, which depends on tissue biomechanical properties, and on the identification of tissue elasticity through the fitting of the experimental displacement images measured using rapid MR-ARFI. Such method was specifically developed for the monitoring of MR-guided High Intensity Focused Ultrasound (MRgHIFU) therapy. Elasticity changes were followed during HIFU ablations (N=6) performed ex vivo in porcine muscle samples, and were compared to temperature changes measured by MR-Thermometry. Shear modulus was found to increase consistently and steadily a few seconds after the heating started, and such changes were found to be irreversible. The shear modulus was found to increase from 1.49±0.48 kPa (before ablation) to 3.69±0.93 kPa (after ablation and cooling). Thanks to its ability to perform quantitative elasticity measurements in a highly localized manner around the focal spot, this method was shown to be particularly interesting for the monitoring of HIFU ablations.

Introduction

Magnetic Resonance-guided High Intensity Focused Ultrasound (MRgHIFU) has emerged over the past decade as a highly valuable non invasive, non ionizing thermal therapy for a large variety of diseases, ranging from oncological [START_REF] Gianfelice | Palliative Treatment of Painful Bone Metastases with MR Imaging-guided[END_REF][START_REF] Schmitz | Image-guided focused ultrasound ablation of breast cancer: current status, challenges, and future directions[END_REF][START_REF] Liberman | Pain Palliation in Patients with Bone Metastases Using MR-Guided Focused Ultrasound Surgery: A Multicenter Study Ann[END_REF] to neurological diseases [START_REF] Zaaroor | Magnetic resonance-guided focused ultrasound thalamotomy for tremor: a report of 30 Parkinson's disease and essential tremor cases[END_REF][START_REF] Elias | A Pilot Study of Focused Ultrasound Thalamotomy for Essential Tremor N[END_REF]. HIFU thermotherapy relies on the absorption of the acoustic energy by the tissue, leading to highly localized temperature elevation. Performing thermal ablations using HIFU therapy requires guidance by medical imaging for the targeting and the monitoring of thermal ablation in real time, as well as for assessing therapy efficiency post-operatively. MRI is particularly appreciated thanks to the constrast in soft tissues it offers intrinsically, to the possibility of imaging several planes freely with identical quality, and thanks to the possibility of monitoring the ablation in real time using MR Thermometry [START_REF] Quesson | Magnetic resonance temperature imaging for guidance of thermotherapy[END_REF].

The monitoring of thermal ablations in real-time is essential in order to ensure that the target tissue is completely ablated while limiting unwanted damage to the surrounding sensitive structures and healthy tissue. MR thermometry based on the temperature dependence of the Proton Resonance Frequency (PRF) (Ishihara et al 1995, Rieke and[START_REF] Rieke | MR thermometry[END_REF] is currently the gold standard in non-fatty soft tissues, since it allows direct visualization of temperature distribution and computation of the accumulated thermal dose (ThD) that serves as the therapy end point when a predefined threshold is reached [START_REF] Sapareto | Thermal dose determination in cancer therapy[END_REF]. However, an uncertainty of a few degrees can lead to an overestimation of the thermal dose by an order of magnitude or more and therefore to an overestimation of the effective ablated zone. In the last 10 years, significant advances have been reported in preclinical research to accelerate and improve spatial coverage of PRF temperature mapping. MR thermometry with an uncertainty of ~1°C on abdominal organs (liver and kidney) has been reported by combining fast acquisition (echo planar imaging with parallel imaging) sequences with online motion compensation and correction of associated magnetic susceptibility artifacts [START_REF] Roujol | Real-time MR-thermometry and dosimetry for interventional guidance on abdominal organs[END_REF][START_REF] Mougenot | Three-dimensional spatial and temporal temperature control with MR thermometry-guided focused ultrasound (MRgHIFU)[END_REF]. Although it is a sophisticated and promising method, PRF-based MR thermometry cannot be performed in fatty tissues. This is particularly problematic for the treatment of breast tumors which are usually surrounded by fat or in the case of liver cancer for patients suffering from adipocytosis [START_REF] Weidensteiner | Stability of real-time MR temperature mapping in healthy and diseased human liver[END_REF]. In these situations, precise monitoring from ThD imaging may become inaccurate. There is thus a clear lack of real-time monitoring methods that could offer additional types of information in complement to temperature/ThD estimate to account for structural changes during ablation. Tissue elasticity has been shown to be a promising biomarker for tissue changes related to heating. Both ultrasound and MR Elastography studies have demonstrated and quantified elasticity changes associated with thermal ablations [START_REF] Mariani | Real time shear waves elastography monitoring of thermal ablation: in vivo evaluation in pig livers[END_REF], Corbin et al 2016[START_REF] Chen | Assessment of in vivo laser ablation using MR elastography with an inertial driver[END_REF][START_REF] Bing | Concurrent ARFI imaging and HIFU ablation using a diagnostic transducer array and ultrasound system with custom beam sequences Ultrasonics Symposium[END_REF][START_REF] Konofagou | Harmonic Motion Imaging (HMI) for Tumor Imaging and Treatment Monitoring[END_REF][START_REF] Wu | Assessment of thermal tissue ablation with MR elastography[END_REF][START_REF] Larrat | MR-guided transcranial brain HIFU in small animal models[END_REF], illustrating the promising value of elasticity imaging methods for the monitoring of thermal ablations.

Magnetic Resonance Elastography (MRE) has been originally developed as a complementary method for the detection of abnormal tissue stiffening [START_REF] Muthupillai | Magnetic resonance imaging of transverse acoustic strain waves[END_REF], based on the fact that most diseases are associated with substantial changes in biomechanical properties. MRE has shown to be highly valuable in a wide spectrum of clinical applications in oncology, hepatology, cardiovascular disease, or neurological diseases [START_REF] Glaser | Review of MR elastography applications and recent developments[END_REF]. Beyond its potential for diagnosis and disease staging, MRE has also been proposed for assessing thermal lesions. Liver tissue stiffening as a result of a laser ablation was demonstrated in vivo by Chen et al. [START_REF] Chen | Assessment of in vivo laser ablation using MR elastography with an inertial driver[END_REF], in an experiment where MRE acquisitions were performed every 2 minutes during laser heating. A similar result was found both in vitro and in vivo by [START_REF] Corbin | Interventional MR elastography for MRI-guided percutaneous procedures[END_REF] on liver tissue using real-time interventional MRE. Both skeletal muscle and turkey breast were shown to be stiffer after a HIFU ablation [START_REF] Wu | Assessment of thermal tissue ablation with MR elastography[END_REF][START_REF] Larrat | MR-guided transcranial brain HIFU in small animal models[END_REF], as opposed to brain tissue which was found to soften after HIFU treatment performed in vivo in rats [START_REF] Larrat | MR-guided transcranial brain HIFU in small animal models[END_REF]. Monitoring HIFU ablations in real-time using MRE would be theoretically possible using realtime methods such as the one proposed by [START_REF] Corbin | Interventional MR elastography for MRI-guided percutaneous procedures[END_REF] for laser ablations. However, such a solution would be particularly challenging and complex in terms of the experimental setup that requires both a mechanical exciter and the HIFU transducer to be positioned around the area of interest.

Magnetic Resonance Acoustic Radiation Force Imaging (MR-ARFI) has been developed initially as a method for spotting the focal zone with high accuracy [START_REF] Mcdannold | Magnetic resonance acoustic radiation force imaging[END_REF]Maier 2008, Larrat et al 2008).

Similarly to MRE, MR-ARFI relies on the encoding of tissue displacement on the MRI phase signal through the use of motion sensitizing gradients (MSG). Tissue displacement is the result of the acoustic radiation force acting at the focal zone. ARFI displacement-based methods are associated with significant limitations for HIFU therapy monitoring: Despite being correlated to tissue biomechanical properties, ARFI displacement cannot be considered as a reliable surrogate for its elasticity. It depends not only on local elasticity, but also on the mechanical properties of the surrounding tissue, on tissue heterogeneity, and on mechanical boundary conditions [START_REF] Vappou | Non-contact, ultrasound-based indentation method for measuring elastic properties of biological tissues using Harmonic Motion Imaging[END_REF]. MR-ARFI methods have been extended to the measurement of mechanical properties: either displacement can be used [START_REF] Liu | Concurrent Visualization of Acoustic Radiation Force Displacement and Shear Wave Propagation with 7[END_REF][START_REF] Dadakova | Optimization of acoustic radiation force imaging: Influence of timing parameters on sensitivity[END_REF][START_REF] Bour | Real-time monitoring of tissue displacement and temperature changes during MR-guided high intensity focused ultrasound[END_REF] or shear wave velocity, which is related to tissue elasticity [START_REF] Liu | Concurrent Visualization of Acoustic Radiation Force Displacement and Shear Wave Propagation with 7[END_REF][START_REF] Souchon | Transient MR elastography (t-MRE) using ultrasound radiation force: Theory, safety, and initial experiments in vitro[END_REF]. Shear wave velocity-based methods involve a spatial averaging effect: a minimal distance (ex: ~10-20 mm [START_REF] Liu | Concurrent Visualization of Acoustic Radiation Force Displacement and Shear Wave Propagation with 7[END_REF]) is needed in order to estimate properly the shear wave velocity. Therefore, these methods are very limited for measuring tissue elasticity specifically at the focal region, which is precisely where tissue structural (and biomechanical) changes are expected to be the most important. A 2015 study by Payne et al. [START_REF] Payne | A simulation technique for 3D MR-guided acoustic radiation force imaging[END_REF] aimed at simulating MR-ARFI displacement profiles using the Green's formalism for a point-source impulse. Computed displacements were found to be in good agreement with experimental ones on phantoms with varying stiffness. However, there was no direct estimation of tissue elasticity and the proposed mechanical information was not used for the monitoring of HIFU therapy. In this study, we propose a new approach for quantitative estimation of tissue elasticity at the focal zone, using the identification of the MR-ARFI profile. Theoretical background of the proposed method is presented and validation of the identification process is performed on numerical phantoms.

Finally, experimental results obtained on ex vivo muscle undergoing HIFU ablations are presented.

Methods

Theoretical background: Relationship between the MR-ARFI profile and tissue mechanical properties

The method used for estimating biomechanical properties is based on a homemade 2D MR-ARFI push simulator that computes the MRI phase profile as a function of experimental parameters. As it has already been described in MR Elastography, a phase shift is observed when the tissue undergoes a displacement denoted as 𝑑 ⃗ simultaneously to a motion sensitizing magnetic field gradient (MSG) 𝐺 ⃗ of the MRI sequence.

In 2D, this can be written in polar coordinates:

𝚽(𝒓, 𝜽) = 𝜸 ∫ 𝑮 ⃗⃗⃗ (𝒕). 𝒅 ⃗ ⃗⃗ (𝒓, 𝜽, 𝒕) 𝑻 𝑴𝑺𝑮 𝟎 𝒅𝒕 (1)
Where TMSG is the duration of the application of the motion-sensitizing gradient, and γ the gyromagnetic ratio of hydrogen. Displacement is assumed to occur only in the direction perpendicular to the image slice.

In practice, this is achieved by choosing a 2D slice orthogonal to the direction of the ARFI push with a MSG applied in the slice direction. In the model, the ARFI push is assumed to be punctual and localized at r = 0. This is obviously an important assumption that will be discussed later in this manuscript, since the actual displacement field is related to the dimensions and shape of the focal spot. In order to take tissue viscoelasticity into account, the displacement resulting from the ARFI push is supposed to follow a first order exponential law:

{ 𝒅(𝒓 = 𝟎, 𝒕) = 𝑨 (𝟏 -𝒆 -𝒕 𝝉 ) 𝒇𝒐𝒓 𝒕 ≤ 𝑻 𝒑𝒖𝒔𝒉 𝒅( 𝒓 = 𝟎, 𝒕) = 𝑨 (𝒆 - 𝒕-𝑻 𝒑𝒖𝒔𝒉 𝝉 ) 𝒇𝒐𝒓 𝒕 ≥ 𝑻 𝒑𝒖𝒔𝒉 (2)
(3)

Where 𝑇 𝑝𝑢𝑠ℎ is the duration of the ARFI push, τ a time constant related to the material's viscoelasticity, and A the maximum displacement. In these equations, the push is assumed to start at t=0. This displacement results in the propagation of a shear wave propagating along the r direction for each θ direction. In other words, the wavefront arrives at a specific position r at a time 𝑡 =

𝑟 𝑐 𝑠

, where cs is the celerity of the shear wave. In a purely linearly elastic medium, the shear wave velocity is directly related to tissue's shear modulus G and tissue density ρ (supposed to be equal to 1000 kg.m -3 in soft tissues) by:

𝒄 𝒔 = √ 𝑮 𝝆 (4)
Therefore, it is possible to compute the displacement field d(r,θ,t) using both the relationships given by eqs.

(2) and ( 3) and the knowledge of tissue elasticity given by its shear wave velocity cs or its shear modulus G. In this study, the tissue is supposed to be locally isotropic, i.e., the shear wave velocity, and therefore the displacement field do not depend on θ. The proposed simulator aims at computing the MRI phase profile ϕ and uses as input the tissue elasticity (given by cs or (G,ρ)), the duration of the ARFI push Tpush , MSG duration TMSG and shape, as well as the size of the spatial window in which the phase profile has to be calculated. Tissue relaxation constant τ is fixed in this study and assumed to be equal to 4 ms. As it will be discussed later in this manuscript, such an assumption is not expected to have a major impact on the elasticity estimate, since the ARFI profiles that are used for the identification are normalized in amplitude.

Although the question of viscoelasticity is highly important from a biomechanical point-of-view, it is beyond the scope of this first study, which is focused on the identification of tissue elasticity only. The maximum push displacement was normalized (A=1). Figure 1 illustrates an example of the phase profile for two different values of cs (cs=1.5 m/s and cs=2 m/s), using bipolar rectangular MSG applied during 10 8 ms, and an ARFI push duration equal to 5 ms. The influence of tissue elasticity can be clearly seen in the phase profile: qualitatively speaking, the size of the ARFI spot increases with tissue stiffness.

Fitting process for elasticity estimation

The fitting process is performed on normalized ARFI profiles centered on the focus, within a window of 51 mm* 51 mm. Experimental ARFI profiles are interpolated by a factor of 2 in order to fit the resolution of the simulated ARFI profiles. The mean quadratic difference between experimental and simulated profiles is used as a cost function for different values of cs. This cost function is minimized using a bounded optimization algorithm (fminbnd, Matlab, Mathworks), using an initial value of 1 m/s. The shear modulus is calculated subsequently using equation (4).

In order to evaluate the identification process, a numerical study was carried out on numerical phantoms with two different levels of added Gaussian noise (SNR = 2, SNR = 5). For each SNR level, tissue stiffness was varied between cs = 1 m/s and cs = 2.4 m/s, which corresponds to a realistic range for shear wave velocity values in soft tissues. All other parameters were set identically for all simulations: ARFI push duration was set to 5 ms, bipolar rectangular MSG were used (total duration = 10 ms), pixel resolution of 1.5 mm, and relaxation constant τ = 4 ms. For each case, 20 numerical phantoms were generated and the shear wave velocity was identified.

Experimental protocol

Ex vivo experiments were performed on a 1.5T preclinical scanner (Siemens Aera, Erlangen, Germany)

combined with an MR-compatible HIFU platform (Image-guided Therapy SA, Pessac, France). This platform consisted of a 256-elements phase array transducer (focal length 14 cm, aperture 13 cm, 1 MHz operating frequency). HIFU sonication was planned and executed from Matlab, using a custom library (IGTfus ,Image-guided Therapy SA, Pessac, France) which controls pulses durations, amplitudes and phases for each element. Two types of experiments were performed:

-The first experimental protocol aimed at validating the proposed MR-ARFI simulator. This was achieved by varying MR-ARFI parameters such as the delay between the MSG and the ARFI push (see TDelay in All experiments were performed on fresh ex vivo pig muscle samples. A Plexiglas tank was positioned on the bed of the HIFU transducer and filled with water at a temperature of 25°C (see Figure 2b). The sample was maintained in a plastic box. Mylar membranes (50 µm thick, smaller than 1/30 of the ultrasonic wavelength) were glued to the base of the plastic box and water tank to ensure ultrasound coupling between the different interfaces. Optimal coil coverage was realized using two 32-element MR coils positioned on each side of the Plexiglas tank and an additional 19 cm loop coil positioned close to the sample.

MR data acquisition was performed using a gradient echo single-shot echo planar imaging (EPI) sequence [START_REF] Bour | Real-time monitoring of tissue displacement and temperature changes during MR-guided high intensity focused ultrasound[END_REF] where a bipolar MSG was added before the EPI echo train (Figure 2a). In the post-treatment stage for tissue elasticity identification, experimental displacement profiles were computed as follows: Phase difference images (φ + -φ -) were unwrapped spatially in 2D using the branch cut method proposed by Goldstein et al. [START_REF] Goldstein R M | Satellite radar interferometry: Two-dimensional phase unwrapping[END_REF], and were normalized with respect to the maximum absolute displacement. The identification process was performed within a limited region around the focal spot (5.1 cm *5.1 cm), on unwrapped, phase difference data. These normalized experimental displacement data were fitted by our MR-ARFI push simulator in order to identify tissue elasticity around the focal zone, as already explained in the numerical phantom study.

Experimental protocol #1: Comparison study

In this first experiment, experimental ARFI profiles were compared to simulated profiles at different values of the delay TDelay between the MSG and the ARFI push. The HIFU shot (SARFI) was synchronized with the MSG with an adjustable delay TDelay set from the beginning of the MSG. The amplitude of the shot (STHERMO)

was turned off to avoid heating (see Figure 2a). The duration and acoustic power of (SARFI) were set to TPush = 2 ms and 447 W, respectively. The delay TDelay was varied between 0 and 9.5 ms by steps of 500 µs, for

Experimental protocol #2: Monitoring of HIFU ablations

In these experiments (N=6) performed on separate muscle samples, a second HIFU shot (STHERMO) was applied during the rest of the sequence, except during the second lobe of the MSG (Fig. 2(a)). An additional adjustable delay TRelax = 20 ms was set before the first lobe of the MSG in order to allow tissue mechanical relaxation after STHERMO and thus avoiding residual displacement encoding.

The sonication was divided in three phases (pre-ablation, ablation, and post-ablation) without interruption of the MR sequence (Fig. 2(c)):

1. Pre-ablation: SARFI = 447 W and STHERMO = 0 W acoustic power, during 30 s. This phase was used to characterize the tissue elasticity properties (young modulus) before ablation with a sonication duty cycle of 1 % resulting in a limited energy deposition.

2. Ablation: SARFI = 447 W and STHERMO = 155 W acoustic power, during 120 s. During this phase, acoustic energy deposition was increased by the use of STHERMO to reach a sonication duty cycle of 86 % in order to induce a thermal lesion.

3. Post-ablation: SARFI = 447 W and STHERMO = 0 W acoustic power, during the rest of the acquisition time. This phase was performed to observe the cooling period and evolution of elasticity after ablation.

Temporal variations of tissue elasticity were compared to those of the temperature measured by MR Thermometry.

Results

Numerical phantoms

Figure 3 illustrates the results for the identification process carried out at two different SNR levels, for shear wave velocity values ranging from cs = 1 m/s and cs = 2.4 m/s. For each value of shear wave velocity, 20 numerical phantoms were generated and the shear wave velocity was estimated. Dots and errorbars denote the average and standard deviation across these 20 cases, respectively.

The identification process was shown to be highly reproducible and reliable, with an average bias less than 1% and an average standard deviation of 4% when SNR = 2, and an average bias less than 1% and an average standard deviation of 2% when SNR = 5.

Experimental results on Muscle

Experimental protocol #1: Comparison between experimental and simulated MR-ARFI profile with varying phase shifts

The delay between the onset of the MSG and the ARFI push was varied with incremental steps of 0.5 ms.

Experimental and simulated normalized MR-ARFI profiles are represented in Figure 4 for delay values of TDelay = 0 ms, 1.5 ms, 3 ms and 4.5 ms. During the procedure, temperature increase measured by MR Thermometry did not exceed 3 °C, resulting in a limited energy deposition. The same profiles can be represented in 1D after averaging along 4 radial directions: These results illustrate the similarity between the experimental and the simulated MR-ARFI profiles for different TDelay, and the capability of the proposed ARFI simulator to properly compute the shape and the dimension of the ARFI pattern. Since the TDelay of zero displays the maximal contrast, such a value was selected for the fitting process during HIFU ablation presented in the rest of the study.

Experimental protocol #2: Monitoring HIFU ablations

The MR-ARFI "spot" was found to increase in size versus time in all of the six HIFU ablations, suggesting that the tissue stiffens as a result of the thermal lesion formation. An example of experimental and corresponding numerical MR-ARFI profiles are represented in Figure 6 at different ablation times. The full dynamic is available as a video in the supplementary material. Temperature elevation during the preablation step was negligible (< 1 °C). Corresponding changes in the elastic modulus and in temperature are represented in Figure 7. It can be clearly seen from these graphs that temperature increases immediately after STHERMO is turned on and that this increase stops once HIFU is turned off. Tissue stiffness behaves differently: stiffening starts later and continues even after ablation is over. This suggests that the formation of the thermal lesion continues even after STHERMO is turned off. In this case, the shear modulus increased from 1.7 kPa (before ablation) to 4.1 kPa after ablation. On average over the six experiments, the shear modulus was found to increase from 1.49±0.48 kPa to 3.69±0.93 kPa after ablation.

Plotting changes in shear modulus against temperature allows observing the different phases of the experiment: First, temperature increases dramatically without any significant changes in stiffness. Stiffness starts increasing slightly towards the end of the ablation. Once ablation is over, temperature starts decreasing, whereas stiffness continues to increase until it reaches its final value.

Discussion

This study proposes a framework that allows computing quantitative MR-ARFI profiles depending on several experimental parameters, such as tissue elasticity, (SARFI) pulse duration, synchronization between the ARFI push and the MSG (TDelay), and properties of MSG.

The fitting process was first evaluated on numerical phantoms with different levels of added noise. Original elasticity was estimated with high accuracy, with low bias and low standard deviation for SNR down to 2.

This numerical study allowed establishing the robustness and the reliability of the identification process.

Simulated MR-ARFI profiles were first compared to experimental ARFI pushes obtained with varying time delays between the ARFI push and MSG. Although the proposed method is effective for any value of TDelay, the fitting process for tissue elasticity identification was found to be the most robust for TDelay = 0 ms, i.e., when the ARFI push starts with the first lobe of the MSG. This may be explained by the fact that the corresponding profile is bipolar (positive and negative peaks, cf Figure 5), which constrains the fitting process and renders it more specific. Very good agreement between experimental and simulated MR-ARFI profiles allowed validating the proposed simulation framework. The MR-ARFI simulator was then used to fit experimental profiles obtained during HIFU ablations on porcine muscle tissue. This allowed identifying changes in tissue elasticity at every pair of MR images received. This study illustrates the value of tissue stiffness as a biomarker for tissue thermal damage. As opposed to temperature, it is directly linked to tissue damage, since tissue stiffness changes remain permanent even after HIFU ablation is turned off. Tissue

stiffness measured with the proposed method appears as a complementary biomarker to thermal dose. The accuracy and the quantitative nature of the proposed method needs to be further investigated, for example through comparison with mechanical testing. However, it is already safe to claim at this stage that it does not suffer from some of the limitations of Thermal Dose, such as its high sensitivity to temperature uncertainty. Comparing these two biomarkers (Thermal Dose and ARFI-derived Elasticity) against histology would be highly interesting and is an obvious continuation of the present study.

Monitoring of thermal ablations requires flow of information in real-time in order to allow immediate action if needed. Although MRI acquisitions were performed fast enough for real-time control of the therapy (~1Hz), all elasticity estimations presented in this study were computed offline. However, translating the identification of the shear modulus in real time is feasible, using an online protocol similar to the one already proposed in interventional MRE [START_REF] Corbin | Interventional MR elastography for MRI-guided percutaneous procedures[END_REF]. The computational time of the identification process varies currently between 1.5 and 3 seconds on a regular computer using the Matlab environment. This identification process should be optimized and shortened in order to allow for its use online in real time.

This study relies on several important assumptions that need to be reminded. First of all, this method assumes that tissue is locally isotropic, which may be questionable in several organs. The computed, simulated MR-ARFI profile does not depend on the direction of the wave propagation, and lead therefore to a circular MR-ARFI profile. Anisotropy would result in a non-circular profile. At this point, our method does not allow identifying tissue elasticity over several directions. However, the proposed MR-ARFI simulator could be extended to anisotropic media in order to identify tissue elasticity against direction of wave propagation G (θ), similarly to what has been proposed in quantitative, anisotropic MRE [START_REF] Chatelin | An automatic differentiation-based gradient method for inversion of the shear wave equation in magnetic resonance elastography: specific application in fibrous soft tissues[END_REF]. The second major assumption behind this study is the nature of the acoustic radiation force that was assumed punctual in our simulator. In reality, the force applies within a focal zone that depends on transducer geometry and properties. In the present study, the width of the focal zone at -3 dB was 1.8 mm in diameter, in the same range of the voxel size of the MR acquisition sequence. Thus, the assumption made on punctual excitation appears reasonable. The higher the quality of focusing is, the closer we are to the assumption of a quasi-punctual force. This assumption may explain the slight differences between the experimental and simulated MR-ARFI profiles (see Figure 4 and5))

In the present study, the monitoring was performed on a single slice orthogonal to the direction of ultrasound propagation, with 5 mm thickness and positioned at the ultrasound focus. Fast single shot EPI sampling allowed us to average acquisitions (five times) to produce phase images of better quality, while maintaining a refresh rate of 1 Hz. Further validation studies aiming at comparing elasticity to thermal dose may require improving the spatial coverage in order to prevent missing out-of-plane heating. From this perspective, already existing methods such as simultaneous multi-slice acquisitions would allow acquiring several slices without changing the proposed acquisition strategy and sequence parameters [START_REF] Barth | Simultaneous multislice (SMS) imaging techniques[END_REF][START_REF] Borman | Towards real-time thermometry using simultaneous multislice MRI[END_REF].

In this study, displacement consecutive to the ARFI push was assumed to follow a 1 st order exponential curve with a relaxation time of 4 ms. This is in the order of magnitude of what has already been measured in previous MR-ARFI studies [START_REF] Dadakova | MR-ARFI for the Quantification of Tissue Elastic Properties[END_REF]. Measuring this relaxation time using the method proposed by [START_REF] Dadakova | MR-ARFI for the Quantification of Tissue Elastic Properties[END_REF] is expected to improve our elasticity estimates, by setting a parameter that is currently unknown. Typical tissue relaxation times are significantly higher under stress relaxation or creep experiments [START_REF] Qiu | Ultrasonic Viscoelasticity Imaging of Nonpalpable Breast Tumors: Preliminary Results[END_REF]. However, in our case, the impulse is particularly short, leading to a particular bandwidth containing high frequencies. Typical relaxation times in the order of the millisecond are not uncommon at high frequency deformation, as attested by the complex modulus spectrum found in soft tissues [START_REF] Papazoglou | Multifrequency inversion in magnetic resonance elastography[END_REF][START_REF] Vappou | Magnetic resonance elastography compared with rotational rheometry for in vitro brain tissue viscoelasticity measurement[END_REF]. In any case, this is a particularly complex issue that is beyond the scope of this study. In the present study, the assumption made on the relaxation time is not expected to influence significantly the identification of tissue elasticity: it will mainly influence the amplitude of the displacement. Since data processing is performed on normalized data, errors on relaxation times are expected to have minimal impact of the results.

The experimental results presented in this study illustrate that elasticity is an interesting biomarker for monitoring HIFU ablations. The general evolution of elasticity is in agreement with previous studies that have investigated changes in tissue elasticity during thermal ablations [START_REF] Arnal | Monitoring of thermal therapy based on shear modulus changes: II. Shear wave imaging of thermal lesions[END_REF][START_REF] Brosses | Temperature dependence of the shear modulus of soft tissues assessed by ultrasound[END_REF][START_REF] Corbin | Interventional Magnetic Resonance Elastography for MRI-guided percutaneous procedures[END_REF]. Whereas tissue temperature was found to increase immediately after HIFU was turned on, elasticity changes seemed to occur after a few seconds. As expected, tissue temperature decreased immediately after HIFU was turned off, as opposed to tissue elasticity, which continued to increase until it reached its final value. This feature of tissue elasticity makes it an interesting biomarker for tissue integrity.

The G vs T curve (Fig. 7) was similar to those reported in the literature [START_REF] Brosses | Temperature dependence of the shear modulus of soft tissues assessed by ultrasound[END_REF][START_REF] Arnal | Monitoring of thermal therapy based on shear modulus changes: II. Shear wave imaging of thermal lesions[END_REF]. As in these studies, very slight softening was found to occur during the initial phase of heating, although this finding should be considered with caution and would require additional experiments to evaluate its consistency. Using tissue biomechanics for monitoring thermal ablations has already been proposed by several research teams. Compared to displacement-based methods (E. [START_REF] Konofagou | Harmonic Motion Imaging (HMI) for Tumor Imaging and Treatment Monitoring[END_REF][START_REF] Bour | Real-time monitoring of tissue displacement and temperature changes during MR-guided high intensity focused ultrasound[END_REF], the proposed method offers the advantage of providing a quantitative parameter directly related to tissue elasticity, as opposed to displacement, which depends highly on many parameters that are not controlled, such as mechanical boundary conditions. The proposed method could be considered as a shear-wave-based elastography method, since it relies of the estimation of the shear wave velocity at the focal region. As opposed to MRE or shear wave-based ultrasonic methods such as shear wave ARFI or Supersonic Shear Wave Imaging (SSI), the proposed method is optimized for local estimation of elasticity specifically at the focal region. This entails significant advantages over the aforementioned methods for specific applications such as monitoring HIFU ablations. However, the proposed method does not allow for mapping tissue elastic properties beyond the focal region, as opposed to MRE. It is expected to perform better at a localized, point-size scale, without the possibility for mapping tissue properties. Beyond its applications for ablation monitoring, the proposed method could also be interesting for quantitative, highly localized biomechanical characterization of soft tissues.

Conclusion

A new, MR-ARFI simulator has been developed that allows computing normalized MR-ARFI profiles.

This tool was used to identify tissue's elastic modulus at the focal region. This identification process was assessed in vitro on muscle tissue and was found to be highly reliable. Elasticity changes were followed during a HIFU ablation, and were compared to temperature changes measured by MR-Thermometry. Shear modulus was found to increase consistently and steadily a few seconds after the heating started, and such changes were found to be irreversible. Thanks to its ability to perform quantitative elasticity measurements in a highly localized manner around the focal spot, this method is particularly interesting for the monitoring of HIFU ablations.

Figure 1 :

 1 Figure 1: (a) Normalized displacement profile used for this example, with an ARFI push duration of 5 ms; (b) MSG profile used (bipolar rectangular); (c) Normalized phase profile obtained with cs = 1.5 m/s; (d) Normalized phase profile obtained with cs = 2 m/s

  Figure 2) and by comparing experimental MR-ARFI profiles to numerical ones. 10 -The second experimental protocol aimed at monitoring HIFU ablations by following the elasticity estimated by the proposed method. Elasticity values were compared to temperature measured by MR-thermometry.

  The acquisition slice was positioned in coronal orientation and MSG direction was set parallel to the acoustic propagation axis to encode longitudinal tissue displacement. In the sequence user interface, duration (TMSG) and amplitude (G, maximal value of 24 mT/m) of MSG were adjustable parameters. In this protocol, the total duration and the amplitude of the MSG cycle were set to TMSG=10 ms and G = 24 mT/m, respectively. This sequence allows simultaneous MR Thermometry and MR-ARFI, similarly to what has already been reported in MRE and MR-ARFI[START_REF] Le | Feasibility of simultaneous temperature and tissue stiffness detection by MRE[END_REF][START_REF] Auboiroux | ARFI-prepared MRgHIFU in liver: Simultaneous mapping of ARFI-displacement and temperature elevation, using a fast GRE-EPI sequence[END_REF].Both experiments followed these steps:1. Localization of muscle sample regarding to the HIFU transducer was performed using a 2D multislice balanced-SSFP sequence: FOV = 250x250 mm², TR/TE/FA = 493 ms/1.36 ms/80 °, 1x1x3 mm 3 voxel size, with a bandwidth of 1149 Hz per pixel. 40 slices of each orientation centered on the transducer were acquired. From this stack of images, a target was selected and adjustments of the position of the monitoring slice were performed.2. One coronal slice with fat saturation was acquired with the proposed sequence. Sequence parameters were: FOV = 189x189 mm², TR/TE/FA = 200 ms/28 ms/35°, voxel size = 1.5x1.5x5 mm 3 , with a bandwidth of 1395 Hz per pixel, GRAPPA acceleration (factor = 2), partial Fourier 6/8, G = 24 mT/m, TMSG = 5 ms for the MSG. Each measurement was averaged five times to increase the precision of phase measurement, resulting in an actual refresh rate of 1 Hz. MR raw data were reconstructed using the open source framework Gadgetron[START_REF] Hansen | Gadgetron: An open source framework for medical image reconstruction[END_REF].The reconstruction process encompasses EPI ghost correction using three acquisition of the central k-space line, coil compression and GRAPPA parallel imaging reconstruction. In both experiments, magnitude and phase images were transferred by TCP/IP to a visualization console (Thermoguide™, Image-guided Therapy SA, Pessac, France) for online computation and display of temperature and displacement maps, using a software module developed in-house[START_REF] Bour | Real-time monitoring of tissue displacement and temperature changes during MR-guided high intensity focused ultrasound[END_REF]. The same console controlled the hardware of the HIFU device for selection of sonication parameters.

  Figure 2: a) Chronogram of the single-shot echo planar imaging MR-ARFI sequence integrating MSG (in blue) with alternating polarities (φ + and φ -). (SARFI) represents the HIFU sonication of duration TPush which produces the ARFI push. TDelay and TMSG represent an adjustable synchronization delay and MSG total duration, respectively. TTL represents the synchronization pulse produced by the scanner. b) 2D balanced-SSFP images acquired in sagittal (top) and transverse (bottom) of the experimental setup. Blue lines represents the acquisition slice centered on the HIFU focal location. c) Typical sonication chronogram combining a pre-ablation, an ablation and a post-ablation phase. Blue box represent successive SARFI shots all along the procedure and additional red box representing STHERMO shots during the ablation phase.

Figure 3 :

 3 Figure 3: Shear wave velocity identified on numerical phantoms with varying elasticity, for two levels of signal to noise ratio. For each shear wave velocity value, 20 numerical phantoms were generated. Dots and errorbars denote the average and standard deviation across these 20 cases, respectively.

Figure 4 :

 4 Figure 4: Experimental and simulated normalized MR-ARFI profiles obtained for 4 different values of TDelay.

Figure 5 :

 5 Figure 5: 1D profiles averaged along radial directions for the data presented in Figure 4, showing the experimental and simulated ARFI profiles for different values of TDelay.

Figure 6 :

 6 Figure 6: Examples of experimental MR-ARFI profiles at different ablation times for case #1, and numerical profiles obtained by fitting the theoretical MR-ARFI model to experimental data. The color scale is the same [-1;1] and corresponds to normalized displacement fields.

Figure 7 :

 7 Figure 7: Example of the evolution of the shear modulus and of the temperature versus time (left); Shear modulus plotted against temperature (right). Timing of the HIFU sonication is represented by the red box

Acknowledgements

This work was partly funded by the French state funds managed by the ANR (within the Investissements d'Avenir programme for the Labex CAMI and TRAIL); Grant number: ANR-11-LABX-0004 and ANR-10-LABX-57, ANR-10-IAHU-04 (IHU Liryc) and by the WorkPackage#3 (Interventional Imaging) of the France Life Imaging (FLI) network.