Kaouther Moussa
email: kaouther.moussa@uphf.fr

Bob Aubouin-Pairault
email: bob.aubouin-pairault@grenoble-inp.fr

Mazen Alamir
email: mazen.alamir@grenoble-inp.fr

Thao Dang

Data-based Extended Moving Horizon Estimation for MISO Anesthesia Dynamics *

This paper presents an extended moving horizon observer to estimate both the states and the pharmacodynamic (PD) parameters of an anesthesia model, based on real data. The inputs of this model are the injection rates of Propofol and Remifentanil. The states represent the concentration of the anesthetic agents in different compartments of the human body (muscles, fat, blood) and in the effect site. The considered output is the Bispectral index (BIS) which is derived from the electroencephalogram (EEG). The observer is designed such that the parameters are estimated during the anesthesia induction phase, and then almost frozen for the rest of the surgery. The estimator is validated on real data that were extracted from the VitalDB database [12].

I. INTRODUCTION

During surgery, the anesthesiologists use different types of anesthetic agents in order to induce and maintain the depth of unconsciousness (hypnosis), the absence of movement (areflexia), and the absence of pain (analgesia) [10].

The level of hypnosis is measured by monitoring the electroencephalogram (EEG) signals and combining several extracted features to determine the depth of anesthesia, usually via processed variables such as the Bispectral Index (BIS) which takes values in the interval [0, 100]. The range [40, 60] characterizes a moderate hypnosis state and the value 50 is generally considered as suitable for surgery.

Compartmental modeling of the drugs distribution has been widely used in the literature, it relies on the conservation principle applied to the exchange of chemicals between different coupled macroscopic systems (compartments) [START_REF] Ionescu | Evaluation of a propofol and remifentanil interaction model for predictive control of anesthesia induction[END_REF]. The anesthetic agent concentrations are assumed to be uniform within a given compartment and the transport rate leaving a compartment is assumed to be proportional to its corresponding concentration.

The 4 th order Pharmacokinetic-Pharmacodynamic (PK-PD) model has been commonly considered in the literature. This model considers four compartments, the central one (blood), the fast equilibrating one (muscles), the slow equilibrating one (fat), and the effect site. In [10] a similar structure has been considered using the Schnider model for Propofol pharmacokinetic [START_REF] Schnider | The Influence of Method of Administration and Covariates on the Pharmacokinetics of Propofol in Adult Volunteers[END_REF] and the Minto model for Remifentanil phamacokinetic [START_REF] Minto | Influence of Age and Gender on the Pharmacokinetics and Pharmacodynamics of Remifentanil: I. Model Development[END_REF].

Estimating the internal states and the parameters of a dynamical model is important to achieve robust performances when using state feedback schemes. For linear unconstrained systems, it is well established that the Kalman filter is the optimal state estimator, when the states and the inputs are subject to normally distributed noises [START_REF] Haseltine | Critical evaluation of extended kalman filtering and moving-horizon estimation[END_REF]. In the presence of hard constraints, such as non-negative drug concentrations, Kalman filtering is not directly applicable. In the case of nonlinear constrained systems, Moving-Horizon-Estimation (MHE) has proven to be a powerful tool for state estimation [START_REF] Schiller | Robust stability of suboptimal moving horizon estimation using an observer-based candidate solution[END_REF]. Furthermore, with the advances of optimization solvers, MHE suits well the context of online estimation based on a window of available past measurements.

The automation of therapy injections has the potential to decrease the clinical workload while maintaining the same care quality [START_REF] Yin | Conflicting interactions in multiple closed-loop controlled critical care treatments: A hemorrhage resuscitation-intravenous propofol sedation case study[END_REF]. Therefore, many closed-loop strategies have been proposed for anesthesia dynamics, for example, [START_REF] Sánchez | Application of model predictive control and moving horizon estimation for the development of closed-loop inhaled anesthetics administration[END_REF] where the authors proposed a Model Predictive Control (MPC) and an MHE for a Single Input Single Output (SISO) model describing the dynamic interactions between the inhaled anesthetic Sevolufrane and the BIS signal. In [START_REF] Da Silva | Online nonlinear identification of the effect of drugs in anaesthesia using a minimal parameterization and bis measurements[END_REF], an extended Kalman filter has been proposed for a minimally parametrized Multiple-Input-Single-Output (MISO) model describing the dynamics between the Propofol and Remifentanil injection rates (inputs) and the BIS as an output, in order to estimate the states and some PK-PD parameters. However, the results rely on approximations and the presented observer is aggressive in the sense that it allows to fit the output with the noise, which means that the parameters change even during the maintenance phase.

Contrary to PK parameters, the PD parameters do not have standard expressions in terms of patient information. According to [11], the uncertainty introduced by the PD parameters shows a more significant influence on the measurable outputs than the one introduced by the PK parameters. Therefore, in this paper, we propose an extended MHE allowing us to estimate both the states and the PD parameters for the complete MISO system with a nonlinear output equation. Furthermore, the cost function is designed such that the penalty on parameters is activated during the maintenance phase in order to avoid having noisy parameters due to the presence of noise on the data, which means that the parameters are identified only during the anesthesia induction phase.

The states being complex to measure practically for the . We show that an online estimation of both the states and the parameters is practically feasible since the average needed time for one iteration varies between 12 and 19ms for a sampling time of 2s. Furthermore, we show that the distributions of the parameters, that were estimated based on real data, closely correspond to practically validated distributions that have been presented in [START_REF] Bouillon | Pharmacodynamic Interaction between Propofol and Remifentanil Regarding Hypnosis, Tolerance of Laryngoscopy, Bispectral Index, and Electroencephalographic Approximate Entropy[END_REF]. This paper is structured as follows: Section II presents the dynamical model that is used for the estimation problem, in Section III the Moving Horizon Estimation technique is recalled and its use for state and parameter estimation is explained. In Section IV we present an observer validation framework. In Section V we present the estimation results for different cases from the VitalDB database and we carry out a statistical analysis on the different estimated parameters. Finally, in Section VI, we present a summary of the results and suggest future perspectives.

II. DYNAMICAL MODEL

Compartmental modeling has been widely used for anesthesia dynamics [10]. The latter consists, firstly, in considering different compartments of the human body, such as the blood, the muscles and the fat, each of them being characterized by its own volume. Then, the concentrations of the anesthetic agents within each compartment are modeled using a set of coupled ordinary differential equations (ODEs).

Pharmacokinetic (PK) allows to characterize the evolution of the injected anesthetic agents in the human body through the absorption phenomenon. Whereas pharmacodynamic (PD) allows to characterize the effect of the drug on the organism, through the BIS in the case of hypnosis study for example.

In the real clinical world, patients usually receive different types of anesthetic agents [START_REF] Yin | Conflicting interactions in multiple closed-loop controlled critical care treatments: A hemorrhage resuscitation-intravenous propofol sedation case study[END_REF]. Therefore, it is important to model the interaction between treatments whether it is synergistic or antagonistic, in order to have a reliable hypnosis indicator, in the case of the BIS for example.

The PK-PD models that are commonly considered for Propofol and Remifentanil are of order 4 (as depicted in Fig. 1) The interaction between both drugs is modeled through the relation between the BIS and the effect site concentration of Propofol (x p4 ) and Remifentanil (x r4 ). We describe in the sequel the different parts of the model that has been considered in this paper.

A. Propofol PK-PD modeling

Propofol is a hypnotic agent that is used to induce and maintain the hypnosis of the patient during surgery. The PK-PD linear model linking the Propofol infusion rate u p to the Propofol concentration in the effect site x p4 is the following:

  ẋp1 ẋp2 ẋp3 ẋp4   =   -a 11p a 12p a 13p 0 a 21p -a 21p 0 0 a 31p 0 -a 31p 0 a 41p 0 0 -a 41p     x p1 x p2 x p3 x p4   +   1 V 1p 0 0 0   u p (1)
where

x p = (x p1 , x p2 , x p3 , x p4 ) T ∈ R 4
+ stands for the vector of Propofol concentrations in the different considered compartments and u p ∈ R + stands for the Propofol infusion rate. By considering that A p is the state matrix and B p = ( 1V 1p , 0, 0, 0) T , Model ( 2) can be condensed in the following equation:

ẋp = A p x p + B p u p (2)
The matrices A p and B p are considered as in equation (2) since the standard model for pharmacokinetic describes the direct effect that the intravenous input u p has on the blood concentration. Furthermore, it models the intercompartmental clearances between the blood compartment, the fat and muscles compartments, as well as the elimination clearance of the drug from the blood compartment. Moreover, the effect compartment is connected to the central compartment (blood) by a first-order rate constant in order to model the delay between a variation in the blood concentration and the beginning of its corresponding effect.

B. Remifentanil PK-PD modeling

Remifentanil is an opioid with a high metabolic clearance, when administered with Propofol, both drugs interact in a synergistic way on hypnosis and analgesia. The PK-PD model of Remifentanil has a similar form as the Propofol one, with different parameter values. Therefore, it can also be described by the following equation:

ẋr = A r x r + B r u r (3)
where x r = (x r1 , x r2 , x r3 , x r4 ) T ∈ R 4 + stands for the vector of Remifentanil concentrations in the different considered compartments and u r ∈ R + stands for the Remifentanil infusion rate.

Since both Propofol and Remifentanil interact separately on the different compartments, the total PK-PD model can be characterized by the following decoupled system:

ẋp ẋr = A p 0 4×4 0 4×4 A r x p x r + B p 0 4×1 0 4×1 B r u p u r (4)
where 0 l×k stands for a zero elements matrix of dimension l × k. By considering that x = (x p , x r ) T ∈ R 8 and u = (u p , u r ) T ∈ R 2 , we can rewrite model (4) as:

ẋ = A T x + B T u (5)
With A T ∈ R 8×8 and B T ∈ R 8×2 being respectively the state and the input matrices.

The parameters of A T and B T are related to the patient information (age, height, weight and gender). Several expressions have been proposed in the literature, for example, Schnider model in [START_REF] Schnider | The Influence of Method of Administration and Covariates on the Pharmacokinetics of Propofol in Adult Volunteers[END_REF] for the Propofol PK parameters (in A p and B p ), Minto model in [START_REF] Minto | Influence of Age and Gender on the Pharmacokinetics and Pharmacodynamics of Remifentanil: I. Model Development[END_REF] for the Remifentanil PK parameters (in A r and B r ) as well as Eleveld model for both Propofol and Remifentanil PK [START_REF] Eleveld | Pharmacokinetic-pharmacodynamic model for propofol for broad application in anaesthesia and sedation[END_REF], [START_REF] Eleveld | An Allometric Model of Remifentanil Pharmacokinetics and Pharmacodynamics[END_REF].

C. Propofol & Remifentanil synergy modeling

The interaction model linking x p4 and x r4 (as depicted in Fig. 1) determines the output equation that characterizes the BIS. The latter is defined as a Hill function depending on an interaction term U as follows:

y = BIS = E 0 -E max U γ 1 +U γ (6)
where E 0 , E max and γ are static parameters representing, respectively, the initial BIS value, the maximal drug effect on the BIS and the slope of the Hill function. The drug interaction term U is defined as follows:

U = U p +U r 1 -β θ + β θ 2 (7)
with:

U p = x p4 c 50p U r = x r4 c 50r θ = U p U p +U r (8)
where c 50p , c 50r and β are also static parameters, standing respectively for the Propofol half-effect concentration, the Remifentanil half-effect concentration and the interaction coefficient.

The PD parameters E max , c 50p , c 50r , γ and β are patient dependent and have to be identified based on the BIS (output) signal. Several works in the literature proposed to simplify this model by considering E max =E 0 and β = 0, for example in [START_REF] Bouillon | Pharmacodynamic Interaction between Propofol and Remifentanil Regarding Hypnosis, Tolerance of Laryngoscopy, Bispectral Index, and Electroencephalographic Approximate Entropy[END_REF]. In order to carry out a consistent comparison with respect to the clinically validated parameters distributions, we consider the same simplification and estimate only the three parameters c 50p , c 50r and γ, we denote by η = c 50p , c 50r , γ T ∈ R 3 the vector of PD parameters to be estimated. Finally, the global system is defined by linear dynamics and a nonlinear output equation as follows:

ẋ = A T x + B T u y = h(x, η) (9)
III. EXTENDED MOVING HORIZON ESTIMATION Extended observers have been presented in the literature [3] and used to estimate exogenous disturbances or the model parameters, in addition to the states. It consists in considering the disturbances (or the parameters) as part of the state vector, with a given dynamical model or invariant dynamics in the case of static parameters.

Since the objective of this paper is to estimate both the states and the PD parameters, the state dynamics can be extended by adding η = 0 to the dynamics by considering x = (x, η) T as follows:

ẋ = ĀT x + BT u y = h( x) (10)
Discretizing model (10) with a sampling period T e , by using the forward Euler method, results in the following extended discrete time system:

xk+1 = A xk + Bu k = Φ( xk , u k ) y k = h( xk ) (11)
where A = I + T e ĀT and B = T e BT , other discretization methods can also be considered. However, since the PK system is relatively slow and in order to induce lower computation time, the Euler method with the associated sampling period of 1 or 2 sec is found appropriate. Moving Horizon Estimation (MHE) is a widely used observation technique [START_REF] Alamir | Partial extended observability certification and optimal design of moving-horizon estimators[END_REF], allowing to recover the states of a given dynamical system, based on the input and the output measurements. In MHE, the current state is estimated by minimizing a cost over a fixed number of past measurements called the estimation horizon and denoted by N MHE . The decision variables being the states over the measurement window, the cost to minimize represents a trade-off between minimizing the output error and ensuring the compatibility of the estimated states with respect to the model and the previous estimated state.

The main advantage of MHE is the explicit handling of the constraints, especially in the context of anesthesia dynamics estimation where the states represent drug concentrations and have to be non-negative. Moreover, the PD parameters that we aim at estimating have also to be non-negative.

The cost function to be minimized at each time step k can be written as follows:

J N MHE (x k | xk-1 , y m , u m ) = k ∑ i=k-N MHE ∥y m i -h( xi )∥ Q + k ∑ i=k-N MHE +1 ∥ xi -Φ( xi-1 , u m i-1 )∥ R(k) (12)
where xk and xk-1 represent, respectively, the state over the estimation horizon (decision variable) and the previous estimated state up to time k -1. The remaining arguments, namely y m and u m represent the output and the input measurements profiles over the estimation horizon, Q and R represent the penalty matrices. Note that the first summation in the cost allows to ensure the fitting with the output while the second one allows to ensure the compatibility of the estimated state with respect to the evolution of the previous estimated one by the dynamical model.

The compatibility cost is penalized with a time dependent weighting diagonal matrix R(k), following a sigmoid function. This allows to identify the parameters during the anesthesia induction phase and to maintain them for the rest of the surgery through a proper penalization. The diagonal elements of R(k) have the following parameterized form:

r j j (k) = θ 1 + θ 2 e -θ 3 e -θ 4 k (13)
where j = 1, . . . , 13. This function allows to have a bounded weight (θ 1 + θ 2 ) for high values of time k, provided that θ 3 , θ 4 > 0. θ 1 is used to penalize the states representing the drugs concentrations within the compartments, whereas θ 2 is used to penalize the parameters during the maintenance phase. For some parameters, a use of θ 1 is required during the induction phase in order to avoid a high variation.

The optimization problem to be solved at each time step k is the following:

min xk J N MHE (x k | xk-1 , y m , u m ) s.t. xi+1 = Φ( xi , u m i ), ∀i = k -N MHE , • • • , k xi ∈ X , ∀i = k -N MHE , • • • , k (14)
Solving ( 14) at each time step k allows to estimate the states (and the parameters) in a receding horizon way, by sliding the measurements window of length N MHE at each time step k. The dynamical constraints have been implemented with the multiple shooting technique allowing to improve the convergence and the sparsity of the optimization problem by lifting it to a higher dimension.

The set X ∈ R 13 allows to define box constraints on the states and the parameters. Regarding the initial guess of the decision variables, a warm start has been implemented, allowing to initialize the states over the estimation horizon by the previous estimated states. At the first time step k = 0, since no estimation is available, the states x (representing the concentrations) are initialized at 0, since u = 0 at the beginning of the surgery, while the parameters η are initialized at an arbitrary pre-tuned PD model.

Since the tuning of the penalty matrices Q and R has a high impact on the estimation quality, in particular in the the context of real data that are corrupted with noise, these matrices are firstly tuned on noisy simulated data, for which we know the real states and parameters, in order to fit them to a realistic scenario and then they are fine-tuned based on real data. The fine-tuning is done in such a way to ensure a smooth BIS profile that do not catch the noise on the output.

IV. OBSERVER VALIDATION

The observability of the non-extended anesthesia system (without parameters) has already been proved in the literature, for example in [START_REF] Nogueira | Positive state observer for the automatic control of the depth of anesthesia-clinical results[END_REF]. However, when extending the system with the parameters dynamics, proving theoretically the observability becomes complex. Therefore, we propose in this section to assess the observability of the extended system based on simulated data for which we know the real states and parameters.

We will consider in the sequel a population of 100 simulated patients, for whom the PD parameters are drawn using the distributions clinically measured in [START_REF] Bouillon | Pharmacodynamic Interaction between Propofol and Remifentanil Regarding Hypnosis, Tolerance of Laryngoscopy, Bispectral Index, and Electroencephalographic Approximate Entropy[END_REF] and presented in Table I, the maximal and minimal bounds are hypothetical values allowing to truncate the normal distribution in order to avoid unrealistic values in the simulations. The PK parameters of Propofol and Remifentanil are computed using Schnider and Minto models [START_REF] Schnider | The Influence of Method of Administration and Covariates on the Pharmacokinetics of Propofol in Adult Volunteers[END_REF], [START_REF] Minto | Influence of Age and Gender on the Pharmacokinetics and Pharmacodynamics of Remifentanil: I. Model Development[END_REF] [START_REF] Bouillon | Pharmacodynamic Interaction between Propofol and Remifentanil Regarding Hypnosis, Tolerance of Laryngoscopy, Bispectral Index, and Electroencephalographic Approximate Entropy[END_REF], * STANDS FOR HYPOTHETICAL VALUES Problem ( 14) has been successively solved, for the 100 simulated cases, by considering N MHE = 20, T e = 1s and using CasADi [START_REF] Andersson | CasADi -A software framework for nonlinear optimization and optimal control[END_REF]. The state constraints can be set based on the solution of ( 2) and ( 3), given the measured inputs u p and u r . Furthermore, several tests on the estimation errors showed that there is no noticeable improvement for estimation horizons that are higher than 20. The computation time for one optimization step is 26 ms, which is reasonable for an online implementation.

The results presented in Figs. 2-3 correspond to one patient and show that in the absence of noise, the output is well fitted and the parameters are well estimated.

Fig. 4 presents the histograms of the mean relative absolute error of the BIS and the states, when no noise affects the BIS, when the BIS is affected with a white noise of variance 1 and when the penalty term R(k) is constant with no noise on the BIS. We can see that in the noise-free case, the estimation errors are centered around 0 and present a low dispersion, which is almost similar for the case where the BIS is corrupted with a white noise, except for the error of x 2 which presents a higher mean error. In the case where R(k) is considered to be constant, we can see that the errors present a higher dispersion even though no noise has been added to the BIS.

Moreover, Fig. 4 shows that the parameters estimation errors are centered around zero and present a low dispersion in the noise-free case, their means remain close to 0 in the noisy case, with a higher dispersion. We can also see that when using a constant penalization R(k), and without any noise, the parameters errors present higher dispersions than in the noisy case.

To summarize, using a simulated based data framework, it is shown that an insight on the observability of the extended system can be gained and the underlying statistics of the estimation errors seems encouraging. The relevance of considering a time-dependent cost in order to freeze the PD parameters, lowering thereby their dispersion is highlighted.

A. Presentation of clinical data

The clinical data used in this paper have been imported from VitaldB [12]. The latter is an open-source database that was created to facilitate the analysis and the monitoring of vital signs during surgery. This database gives access to highresolution multi-parameter data from 6388 cases with several indicator signals.

We considered a population of 50 cases from the database, this population has a median age of 55, a median height of 64, a median weight of 166 and a 32% proportion of women. They have been selected based on two criteria, the first one is the use of both Propofol and Remifentanil and the second one is the BIS stability. Indeed, the database contains some cases with an unstable BIS during the maintenance phase, which is due to different possible exogenous disturbances. Therefore, in order to avoid skewing the estimation of the parameters with exogenous signals, we consider only the cases with a stable BIS during the maintenance phase. Furthermore, the data have been processed in order to remove the NaNs.

B. MHE implementation on Real-data

The MHE optimization problem has been solved using CasADi [START_REF] Andersson | CasADi -A software framework for nonlinear optimization and optimal control[END_REF] with N MHE = 20 and T e = 2s. The initial guesses were set to the nominal values. The PK parameters have been considered according to the Eleveld Model [START_REF] Eleveld | An Allometric Model of Remifentanil Pharmacokinetics and Pharmacodynamics[END_REF], [START_REF] Eleveld | Pharmacokinetic-pharmacodynamic model for propofol for broad application in anaesthesia and sedation[END_REF], since it showed a better fitting of clinically measured PD parameters distributions presented in [START_REF] Bouillon | Pharmacodynamic Interaction between Propofol and Remifentanil Regarding Hypnosis, Tolerance of Laryngoscopy, Bispectral Index, and Electroencephalographic Approximate Entropy[END_REF]. E 0 has been identified from the first samples of the BIS, since it represents the BIS value before the induction phase. Figs. 56show the results of the extended MHE for the case ID 167, we can see that the parameters are quasi static during the maintenance phase and the estimated BIS fits the measured one even though the latter presents a relatively high variance noise. We can also see that the BIS that was calculated for nominal PD parameters presents a high fitting error, which highlights the importance of estimating the PD parameters. Fig. 7 shows the distribution of the estimated parameters with the practically validated ones, we can see that the distribution of c 50p presents a good fitting, while the distributions of c 50r and γ are slightly shifted. This can be explained by the difference of the population that have been considered for the study in [START_REF] Bouillon | Pharmacodynamic Interaction between Propofol and Remifentanil Regarding Hypnosis, Tolerance of Laryngoscopy, Bispectral Index, and Electroencephalographic Approximate Entropy[END_REF], since the median age presents a considerable difference.

VI. CONCLUSION

We presented in this paper a data-based extended moving horizon estimator for anesthesia dynamics, the latter allows to estimate both the state and the parameters. The estimated variables can be used in a state feedback scheme in order to stabilize the BIS in the interval [40,60]. Moreover, a global statistical analysis can be carried out based on the estimated parameters in order to study the correlation between the different variables. Furthermore, fault detection based observers can be considered in order to detect the presence of external disturbances on the BIS signal. Histograms of the estimated parameters with the practically validated distribution in [START_REF] Bouillon | Pharmacodynamic Interaction between Propofol and Remifentanil Regarding Hypnosis, Tolerance of Laryngoscopy, Bispectral Index, and Electroencephalographic Approximate Entropy[END_REF] surgical patients. Sci Data, 9(1):279, Jun 2022.

Fig. 1 .

 1 Fig. 1. Compartmental scheme of the pharmacokinetic and the pharmacodynamic model of Propofol and Remifentanilcase of anesthesia, the observer is first validated on simulated data, for which the state profiles and the parameter values are available. Then, the observer is applied to different real clinical cases from the VitalDB database[12]. We show that an online estimation of both the states and the parameters is practically feasible since the average needed time for one iteration varies between 12 and 19ms for a sampling time of 2s. Furthermore, we show that the distributions of the parameters, that were estimated based on real data, closely correspond to practically validated distributions that have been presented in[START_REF] Bouillon | Pharmacodynamic Interaction between Propofol and Remifentanil Regarding Hypnosis, Tolerance of Laryngoscopy, Bispectral Index, and Electroencephalographic Approximate Entropy[END_REF].This paper is structured as follows: Section II presents the dynamical model that is used for the estimation problem, in Section III the Moving Horizon Estimation technique is recalled and its use for state and parameter estimation is explained. In Section IV we present an observer validation framework. In Section V we present the estimation results for different cases from the VitalDB database and we carry out a statistical analysis on the different estimated parameters. Finally, in Section VI, we present a summary of the results and suggest future perspectives.

Fig. 2 .Fig. 3 .

 23 Fig. 2. The drug injection rate profiles, the measured and the estimated BIS for N MHE = 20

Fig. 4 .

 4 Fig. 4. Histograms of the Mean Relative Absolute Errors (MRAEs) of the BIS and the states (up) and histograms of the true parameters and their corresponding estimation errors ε ĉ50p , ε ( ĉ50r ) and ε ( γ) (bottom).

Fig. 5 .

 5 Fig. 5. Case ID 167: the drug injection rates (up), the measured and the estimated BIS (down)

Fig. 6 .

 6 Fig. 6. Case ID 167: the estimated parameters

  . .

	Parameter	c 50p	c 50r	γ
	Min *	2	6	1
	Max *	9	25	5
	Mean	4.47	19.3	1.43
	Standard deviation 1.34 5.79 0.73

TABLE I THE

 I PD PARAMETERS STATISTICS CONSIDERED FOR THE OPTIMIZATION PROBLEM

Univ. Grenoble Alpes, CNRS, Grenoble INP, VERIMAG, 38000 Grenoble, France thao.dang@univ-grenoble-alpes.fr

* This work has been partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) funded by the french program Investissement d'avenir, the RITMEA research program funded by the European Regional Development Fund, the State, and the Hauts-de-France Region as well as the Regional Delegation for Research and Technology, the Ministry of Higher Education and Research and the CNRS.


Karine Choquet

Jonathan Vappou

Paolo Cabras

Ounay Ishak

Afshin Gangi

Élodie Breton

Elodie Breton

Magnetic Resonance Acoustic Radiation Force Imaging (MR-ARFI) for the monitoring of High Intensity Focused Ultrasound (HIFU) ablation in anisotropic tissue

Keywords: Magnetic Resonance Guided Interventional Procedures, High-Intensity Focused Ultrasound Ablation, Acoustic Radiation Force Impulse Imaging, Anisotropy, Tissue Elasticity Imaging

Introduction

High Intensity Focused Ultrasound (HIFU) thermal therapy is a non-ionizing and non-invasive clinical procedure that relies on localized absorption of the acoustic energy. To ensure the safety and efficiency of the therapy, it is usually performed either under magnetic resonance (MR) [START_REF] Jolesz | MRI-Guided Focused Ultrasound Surgery[END_REF][START_REF] Kim | Advances in MR image-guided high-intensity focused ultrasound therapy[END_REF][START_REF] Kuroda | MR techniques for guiding high-intensity focused ultrasound (HIFU) treatments[END_REF] or diagnostic ultrasound (US) imaging guidance [START_REF] Ebbini | Ultrasound-guided therapeutic focused ultrasound: Current status and future directions[END_REF]. Several non-invasive MR and US temperature monitoring methods have been proposed for estimating the temperature variation during HIFU ablation [START_REF] Ebbini | Real-time Ultrasound Thermography and Thermometry[END_REF][START_REF] Geoghegan | Methods of monitoring thermal ablation of soft tissue tumorsa comprehensive review[END_REF][START_REF] Lewis | Thermometry and ablation monitoring with ultrasound[END_REF][START_REF] Odéen | Magnetic resonance thermometry and its biological applications -Physical principles and practical considerations[END_REF][START_REF] Rieke | MR Thermometry[END_REF]. Among them, the Proton Resonance Frequency Shift (PRFS) MR technique [START_REF] Poorter | Noninvasive MRI thermometry with the proton resonance frequency method: study of susceptibility effects[END_REF][START_REF] Ishihara | A precise and fast temperature mapping using water proton chemical shift[END_REF] is currently the clinical gold-standard [START_REF] Kokuryo | Recent technological advancements in thermometry[END_REF] due to its tissue-independency over a wide temperature range. However, the limitations of PRFS MR-thermometrysuch as its sensitivity to magnetic field variations and the presence of fat moleculescall for the development of complementary biomarkers.

In addition to temperature information, the alteration of mechanical properties related to actual thermal tissue damage have been investigated using elastography techniques [START_REF] Arnal | Monitoring of thermal therapy based on shear modulus changes: II. Shear wave imaging of thermal lesions[END_REF][START_REF] Chen | Assessment of in vivo laser ablation using MR elastography with an inertial driver[END_REF][START_REF] Corbin | Interventional MR elastography for MRI-guided percutaneous procedures: Interventional MR Elastography[END_REF][START_REF] Mariani | Real time shear waves elastography monitoring of thermal ablation: in vivo evaluation in pig livers[END_REF][START_REF] Brosses | The link between tissue elasticity and thermal dose in vivo[END_REF]. Several MR Imaging (MRI)-based methods have been proposed to track biomechanical changes in tissues undergoing HIFU therapy. Using harmonic MR Elastography (MRE) [START_REF] Muthupillai | Magnetic resonance elastography by direct visualization of propagating acoustic strain waves[END_REF], local irreversible increase in stiffness was reported as a result of HIFU ablation in ex vivo porcine muscle tissues [START_REF] Wu | Assessment of thermal tissue ablation with MR elastography[END_REF]. It is worth noting that, PRFS temperature information may also be recovered simultaneously from MRE data using a gradient echo MR sequence [START_REF] Corbin | Interventional MR elastography for MRI-guided percutaneous procedures: Interventional MR Elastography[END_REF][START_REF] Kim | Simultaneous fat-referenced proton resonance frequency shift thermometry and MR elastography for the monitoring of thermal ablations[END_REF][START_REF] Le | Feasibility of simultaneous temperature and tissue stiffness detection by MRE[END_REF].

Using real-time MRE and MR-thermometry, Kim et al. [START_REF] Kim | Simultaneous fat-referenced proton resonance frequency shift thermometry and MR elastography for the monitoring of thermal ablations[END_REF] reported a significant decrease in the shear modulus of ex vivo porcine muscle tissue as a consequence of moderate temperature elevation induced by HIFU. Nonetheless, MRE requires induction of shear waves with an external mechanical exciter that complexifies the interventional set-up, which is usually already bulky. For this reason, using the HIFU transducer as the source of both treatment and shear wave elastography waves represents a promising solution.

Due to absorption and scattering effects in biological tissues, HIFU pulses exert an acoustic radiation force, i.e. tissues locally undergo a micrometric displacement that causes the propagation of shear waves [START_REF] Chu | Acoustic radiation pressure produced by a beam of sound[END_REF][START_REF] Torr | The acoustic radiation force[END_REF]. Derived from the ultrasound Acoustic Radiation Force Impulse (ARFI) Imaging [START_REF] Nightingale | Shear-wave generation using acoustic radiation force: in vivo and ex vivo results[END_REF], the MR-ARFI imaging approach was proposed [START_REF] Larrat | Radiation force localization of HIFU therapeutic beams coupled with magnetic resonance-elastography treatment monitoring in vivo application to the rat brain[END_REF][START_REF] Mcdannold | Magnetic resonance acoustic radiation force imaging[END_REF], based on encoding the acoustic displacement in the MR phase signal using an appropriate motionsensitive sequence. Originally, MR-ARFI consisted in generating a brief ultrasound wave train and then imaging this displacement to locate the focal spot or to correct for acoustic aberrations for HIFU therapy planning [START_REF] Holbrook | In vivo MR acoustic radiation force imaging in the porcine liver[END_REF][START_REF] Kaye | Rapid MR-ARFI method for focal spot localization during focused ultrasound therapy[END_REF][START_REF] Larrat | MR-guided transcranial brain HIFU in small animal models[END_REF]. The use of ARFI has also been suggested to perform MR transient elastography [START_REF] Souchon | Transient MR elastography (t-MRE) using ultrasound radiation force: Theory, safety, and initial experiments in vitro[END_REF], paving the way for monitoring HIFU ablations. As ARFI-induced tissue displacements are inversely correlated to the stiffness of the tissue, multiple papers described MR-ARFI displacement-based monitoring methods [START_REF] Auboiroux | ARFI-prepared MRgHIFU in liver: simultaneous mapping of ARFI-displacement and temperature elevation, using a fast GRE-EPI sequence[END_REF][START_REF] Bever | Simultaneous MR thermometry and acoustic radiation force imaging using interleaved acquisition[END_REF][START_REF] Bitton | Toward MR-guided high intensity focused ultrasound for presurgical localization: Focused ultrasound lesions in cadaveric breast tissue[END_REF][START_REF] Bour | Real-time monitoring of tissue displacement and temperature changes during MR-guided high intensity focused ultrasound[END_REF][START_REF] Ilovitsh | Acoustic radiation force imaging using a single-shot spiral readout[END_REF][START_REF] Qiao | Simultaneous acoustic radiation force imaging and MR thermometry based on a coherent echo-shifted sequence[END_REF]. However, ARFIinduced displacements not only depend on the local tissue elasticity, but also on mechanical boundary conditions and tissue acoustic absorption, the latter parameters being susceptible to change during thermal ablations. To overcome this limitation, Vappou et al. [START_REF] Vappou | MR-ARFI-based method for the quantitative measurement of tissue elasticity: application for monitoring HIFU therapy[END_REF] developed a method that relates MR-ARFI spot profile to tissue's local shear wave velocity and, by consequence, assuming a linear viscoelastic model, shear modulus during HIFU ablations. The paper reported an irreversible increase in tissue stiffness of ex vivo HIFU-ablated tissue, making it an interesting biomarker of thermal damage. However, the linear viscoelastic model relied on two strong hypotheses: 1-the viscoelastic relaxation time  was arbitrarily set from the literature, hence not tissue-specific, 2-the tissue was assumed to be isotropic, which may be inaccurate for several soft tissues such as muscle or kidney. Indeed, in such anisotropic media, shear waves generated for elastography purposes are likely to propagate with different velocities depending on their direction [START_REF] Schmidt | Magnetic resonance elastography of slow and fast shear waves illuminates differences in shear and tensile moduli in anisotropic tissue[END_REF].

In this study, we propose an MR-ARFI-based method that provides both temperature maps and stiffness polar maps for the monitoring of HIFU therapy. The viscoelastic model takes the 5/27 tissue-specific  into account, and the geometric attenuation of the mechanical wave. A validation experiment is performed in a homogeneous phantom calibrated for elastography. The capability of the proposed approach to monitor a HIFU ablation is tested in a tissue-mimicking gel phantom as well as in ex vivo and in vivo swine muscle tissue.

Materials and Methods

MR-ARFI for HIFU ablation monitoring

Following the principle of MR-ARFI first proposed by McDannold and Maier [START_REF] Mcdannold | Magnetic resonance acoustic radiation force imaging[END_REF], the HIFU transducer delivering the HIFU therapy is also used to generate the ARFI mechanical wave within the tissue to be monitored. For that purpose, a short sonication pulse of the order of one millisecond is applied. This generates an acoustic radiation force which induces a micrometric displacement (push) in the focal zone (Fig. 1a), resulting in the propagation of a transient shear wave of velocity cs. This mechanical wave is encoded on the MRI phase signal using bipolar motion sensitizing gradients (MSG) synchronized to the acoustic radiation force (Fig. 1b). In general, stiffer tissues lead to a wider ARFI spot.

Vappou et al. [START_REF] Vappou | MR-ARFI-based method for the quantitative measurement of tissue elasticity: application for monitoring HIFU therapy[END_REF] proposed to identify the experimental MR-ARFI phase profile to a theoretical one, allowing determination of cs and therefore the local shear modulus μ at the focal spot, assuming that the medium is isotropic, linear viscoelastic and locally homogeneous of density 𝜚  1000 kg / m 3 :

𝜇 = 𝜚𝑐 𝑠 2 Equation 1
The theoretical MR-ARFI profile in the phase image is given by:

𝜙 𝑡ℎ (𝑟,  ) = 𝛾 ∫ 𝐺(𝑡)𝛿(𝑟,  , 𝑡)𝑑𝑡 𝑇 𝑀𝑆𝐺 0 Equation 2
where  is the gyromagnetic ratio of the hydrogen nucleus, G is the MSG amplitude, TMSG its duration, and δ(r,,t) is the displacement at time t written in polar coordinates (r, ) (Fig. 2). To model tissue viscoelasticity, δ at the focal spot (r = 0) is assumed to follow a 1 st order exponential law of time constant  representing the mechanical relaxation time of the tissue:

𝜹(𝒓 = 𝟎, 𝜽, 𝒕) = { 𝟎 ≤ 𝒕 ≤ Tpush : 𝜹(𝒕) = 𝑫 (𝟏 -𝒆 -𝒕 𝝉 ) t ≥ Tpush : 𝜹(𝒕) = 𝑫 (𝟏 -𝒆 - 𝑻 𝒑𝒖𝒔𝒉 𝝉 ) 𝒆 - 𝒕 -𝑻 𝒑𝒖𝒔𝒉 𝝉 Equation 3
where τ is the mechanical relaxation time of the tissue, Tpush corresponds to the duration of the ARFI push and D to the asymptotic value of the ARFI displacement at the focal spot. The ARFI displacement induced away from the focal spot is described along radial propagation lines (r > 0 at a given  angle) considering the medium as locally homogeneous:

𝛿(𝑟, 𝜃, 𝑡) = { 𝑡 < 𝑟 𝑐 𝑠 : 𝛿(𝑡) = 0 t ≥ 𝑟 𝑐 𝑠 : 𝛿(𝑡) = 1 √𝑟 𝛿(0, 𝜃, 𝑡 - 𝑟 𝑐 𝑠 ) Equation 4
with cs denoting the shear wave velocity. Hence, the propagation of the shear wave is modeled as a rectilinear translation along 1D radial lines of the displacement δ( r = 0, , t) induced at the focal point (Fig. 2), where r c s is the time delay corresponding to the arrival of the wavefront at position r. Assuming a cylindrically-shaped wavefront, the attenuation in amplitude is modeled as inversely proportional to the square root of the distance (r) from the excitation source.

Anisotropic MR-ARFI for HIFU ablation monitoring

In this paper, the approach is further developed to estimate tissue anisotropy around the focal spot. The shear modulus μ is computed in predetermined radial directions around the focal spot, by individually fitting 1D radial phase profiles to the linear viscoelastic model assuming that the medium is locally homogeneous along each radial direction (Fig. 2).

9/27

The specific relaxation time  of the tissue at the focal spot is determined experimentally prior to the HIFU treatment following the approach published by Kaye et al. [START_REF] Kaye | Adapting MRI Acoustic Radiation Force Imaging For In Vivo Human Brain Focused Ultrasound Applications[END_REF] and Dadakova et al. [START_REF] Dadakova | Optimization of acoustic radiation force imaging: Influence of timing parameters on sensitivity[END_REF]. For that purpose, several MR-ARFI images are acquired with different time shift toffset values (Fig. 1b), and the MR phase signal measured in the central pixel of the focal spot is fitted (Fig. 2, bottom right) according to Equation 2 and 3.

The experimental MR-ARFI images were resampled using linear spatial interpolation by a factor of two for the identification process. Background voxels are first masked out based on their magnitude. Hence, the radial 1D phase profiles centered on the focal spot are extracted from the masked spatially interpolated MR-ARFI phase image. A bounded optimization algorithm [START_REF] Errico | Matlab functions: fminsearchbnd, fminsearchcon[END_REF] is used to maximize the correlation between the experimental and simulated 1D MR phase profiles. The solution domain is constrained according to realistic physical limits, namely 0.8 kPa ≤ μ ≤ 100 kPa for the shear modulus [START_REF] Wu | Assessment of thermal tissue ablation with MR elastography[END_REF][START_REF] Farrer | Characterization and evaluation of tissue-mimicking gelatin phantoms for use with MRgFUS[END_REF].

MR-ARFI experiments

Two different series of experiments were planned. The first experiments aimed at evaluating the accuracy and quantitative nature of the proposed method for elastography on a calibrated phantom. Then, the ability of the proposed method to monitor HIFU heating in different scenarios (phantom, ex-vivo and in-vivo) was evaluated.

The same MR pulse sequence was used for all experiments: a segmented echo planar imaging (EPI) gradient echo sequence with MSG and focused ultrasound (FUS) trigger adjuncts was developed. The toffset parameter is used to trigger the ARFI push relatively to the beginning of the MSG (Fig. 1b). Experiments were performed on 1. thermometry by averaging the two phase images [START_REF] Corbin | Interventional MR elastography for MRI-guided percutaneous procedures: Interventional MR Elastography[END_REF][START_REF] Le | Feasibility of simultaneous temperature and tissue stiffness detection by MRE[END_REF]. PRFS with B0-drift correction is used to measure temperature changes. Phase difference images were spatially unwrapped for elastography, while average phase images were temporally unwrapped for thermometry.

Assuming the shear wave propagates perpendicularly to the US beam axis, the MSG direction was set parallel to this axis to encode axial tissue displacement, and the acquisition slice was positioned perpendicular to this axis (Fig. 1a). A 1 MHz 128-element spherical transducer (Imasonic, France) with a nominal focal length of 6 cm was used for all the experiments, for both ARFI push and therapy. The transducer was driven by an MR-compatible, 256-channel HIFU generator (Image Guided Therapy, Pessac, France). A latex balloon filled with degassed water was used for acoustic coupling. A granular jamming positioning system held the transducer in place during all experiments [START_REF] Cabras | A new versatile MR-guided highintensity focused ultrasound (HIFU) device for the treatment of musculoskeletal tumors[END_REF]. The ARFI push duration was set equal to TMSG/2 = 5 ms.

The quantitative nature of the proposed method was first evaluated in a homogeneous commercial phantom calibrated for elastography (CIRS, Norfolk, USA). The value of reference provided by the manufacturer for shear elastic modulus is 5.9 kPa ± 5%. The acoustic power of the ARFI shots was set to 200 W. The delay toffset between the ARFI push and the MSG was varied from -6 to 1 ms, every 1 ms, these measures being repeated 4 times. These data were used to estimate the viscoelastic relaxation time  at the focus, and to study the sensitivity of the local radial shear modulus μ to toffset.

HIFU heating was then monitored using the proposed method. It was first evaluated in a tissuemimicking phantom, made of porcine gelatin powder (160 g), water (2.6 kg) and condensed milk (800 g) containing 8 % fat and 55% sugar in order to increase the acoustic absorption [START_REF] Farrer | Characterization and evaluation of tissue-mimicking gelatin phantoms for use with MRgFUS[END_REF].

Before refrigeration, the mixture was degassed using a vacuum chamber.

For the ex-vivo experiment, a porcine muscle was obtained from local butcher, whereas the invivo tests were performed on one swine model (male, 32 kg). The pig was anesthetized (1.8% isoflurane) and placed supine with the HIFU transducer positioned on its right thigh.

The FUS sequence consisted in interleaving a HIFU heating phase with an ARFI shooting phase within each TR period (Fig. 3).

Fig. 3 Prior and after each HIFU experiment, the phantom/tissue mechanical relaxation time τ was determined with toffset varying from -5 ms to 4 ms every 1 ms, -2.5 ms to 2 ms every 0.5 ms, and -5 ms to 0.5 ms every 0.5 ms, for the gelatin phantom, ex vivo and in vivo muscle experiments, respectively.

Results

Calibrated phantom experiment

The value of reference determined for  is 3.4 ms at the focal spot. MR-ARFI phase profiles acquired with 4 increasing values of toffset can be seen from left to right in Fig. 4 (top row).

Corresponding μ polar maps averaged over the 4 repetitions are displayed on the bottom row.

13/27

All phase profiles and μ polar maps can be found in a separate figure (Online Resource 1). The obtained polar representations of μ appear circular, characteristic of an isotropic medium. Assuming strict isotropy, the average and standard deviation of µ across all directions  and 4 repetitions for each toffset are given in Table 2. The average across all directions , repetitions and toffset is equal to 8.2 kPa ± 7 % for μ.

MR-ARFI monitoring experiments

For each experiment, the acquisition time for both temperature and μ polar maps was 4 s, corresponding to the acquisition of 2 MR images with opposite motion sensitivity polarities.

Tissue-mimicking phantom

During the HIFU ablation in the gelatin phantom, the MR-ARFI spot diameter was observed to decrease versus time as shown in the video (Online Resource 2), suggesting that the gel softens as a result of heating. This can be observed from the μ polar maps obtained preand postablation (Fig. 5), their circular shape reflecting the high isotropy of the phantom. The shear modulus μ was found to decrease from around 5 kPa to 3 kPa, while temperature increased by 12°C at the central pixel of the focal spot (Fig. 5). Regarding , its value was found constant pre-and post-HIFU at 5.2 ms.

Ex vivo tissue

In the ex vivo porcine tissue experiment, the muscle tissue appeared to be highly anisotropic, as depicted by the oval MR-ARFI spots and the resulting peanut-shaped μ polar maps shown in Fig. 6 before and after HIFU.  increased slightly from 2.2 ms pre-heating to 2.6 ms postheating, while μ globally decreased during HIFU heating (video Online Resource 3). In this HIFU experiment, the temperature increased by 22°C at the central pixel of the focal spot.

In vivo tissue

In porcine muscle in vivo, the MR-ARFI spot was observed to widen progressively with time (video Online Resource 4), which is consistent with the increase in shear modulus μ(θ) observed in the polar plots during the ablation.  was found to increase as well, from 5.8 ms pre-ablation to 13.9 ms post-ablation. The maximum temperature rise reached at the focus was 50 ºC.

Looking at the MR-ARFI phase profiles, one can notice that the images are obscured by an oblique band of phase wrapping-related artifacts. As it was poorly corrected by unwrapping, we excluded this area from the identification process. Fig. 7 then shows the shear modulus angular distribution estimated in the left half of the MR-ARFI spot, both pre-and post-heating in tissues that appeared to be anisotropic.

Discussion

In this paper, we introduce a method that provides both shear modulus polar maps and temperature maps for the monitoring of HIFU therapy. The present MR-ARFI shear wave velocity-based approach allows for the characterization of tissue anisotropy while taking into account the viscoelastic relaxation time  determined before the ablation. Concurrent monitoring and therapy were performed by interleaving short ARFI shots and HIFU heating at the focus during the HIFU treatment. The direction-dependant shear modulus μ is then identified for a user-determined number of 1D MR-ARFI phase profiles centered on the HIFU focus, the result of the identification process being represented in a shear modulus polar map updated during the treatment.

In the commercial elasticity phantom, good repeatability was found between the 4 repetitions and all the polar maps radii of MR-ARFI for each individual toffset. The profile of the shear modulus polar map was consistent with the expected isotropy of the phantom. Although the shear modulus determined with MR-ARFI (8.2 kPa  7 %) was in the same range as the one reported by the manufacturer (5.9 kPa  5 %), there was a significant difference that may be due to differences between mechanical testing protocols. Indeed, the value provided by the manufacturer has been obtained through mechanical indentation whereas our approach is based on the propagation of a transient shear wave. It is well known that quasi-static and dynamic elastography methods may yield different results due to the differences in excitation frequencies [START_REF] Lefebvre | Comparison of Viscoelastic Property Characterization of Plastisol Phantoms with Magnetic Resonance Elastography and High-Frequency Rheometry[END_REF][START_REF] Vappou | Magnetic resonance elastography compared with rotational rheometry for in vitro brain tissue viscoelasticity measurement[END_REF].

Our results highlight the influence of the toffset parameter on the shear modulus identification process. This parameter allows for controlling the size of the MR-ARFI pattern for a given Tpush/TMSG configuration. On the one hand, a small ARFI spot -compared to the spatial resolution of the image -can lead to errors in the fitting process due to the limited amount of 21/27 data available in the radial phase profile. On the other hand, a large ARFI spot is more sensitive to the hypothesis of local homogeneity of the shear modulus along each propagation radius. In addition, a large focal spot may be particularly sensitive to shear wave dispersion effects that are not taken into account in the current model, as the displacement profile is assumed to be unchanged in shape during propagation (Fig. 2 and Eq.4). These extreme values of toffset may result in significant errors in shear modulus estimates, as a result of inaccuracy of the MR-ARFI phase identification process. For these reasons, an optimal value for toffset should be carefully chosen before each monitoring experiment, in order to control the size of the initial MR-ARFI spot.

Resulting from the identification process of each individual radial phase profile of the ARFI spot, the polar maps obtained with MR-ARFI reflected the isotropy of gel phantoms, and the anisotropy of biological soft tissue as observed in ex vivo and in vivo porcine muscle. It should be noted that the number and values of  angles shown in the polar map can be freely defined, taking the spatial resolution and the computation time into account. During HIFU treatment, the polar maps obtained with MR-ARFI mainly reflect the intrinsic anisotropy of the tissue.

They may also help in detecting direction-dependent changes in tissue elasticity during the HIFU treatment, as a result of vascular perfusion for instance. The current model proposes an increment towards the characterization of anisotropic tissues; however, it relies on the strong hypothesis that the medium is radially homogeneous within the ARFI-spot. This hypothesis may be particularly inaccurate in tissues featuring curved muscle fibers, such as the heart.

The decrease of the shear modulus observed during HIFU in the gelatin phantom reflects the gel softening due to HIFU heating. Ex vivo porcine muscle tissue was observed to soften, as already reported in the literature [START_REF] Wu | Assessment of thermal tissue ablation with MR elastography[END_REF][START_REF] Kim | Simultaneous fat-referenced proton resonance frequency shift thermometry and MR elastography for the monitoring of thermal ablations[END_REF][START_REF] Kruse | Tissue characterization using magnetic resonance elastography: preliminary results\ast[END_REF] with moderate temperature increase. The in vivo HIFU ablation in porcine thigh muscle resulted in local tissue stiffening as expected for higher temperature elevations [START_REF] Arnal | Monitoring of thermal therapy based on shear modulus changes: II. Shear wave imaging of thermal lesions[END_REF][START_REF] Wu | Assessment of thermal tissue ablation with MR elastography[END_REF][START_REF] Vappou | MR-ARFI-based method for the quantitative measurement of tissue elasticity: application for monitoring HIFU therapy[END_REF][START_REF] Brosses | Temperature dependence of the shear modulus of soft tissues assessed by ultrasound[END_REF].

Other MR-ARFI-based elastography approaches have been recently introduced. With the combination of multiple ARFI pushes in a plane perpendicular to the ultrasound propagation axis, Hofstetter et al.'s time of flight-based approach [START_REF] Hofstetter | Efficient shear wave elastography using transient acoustic radiation force excitations and MR displacement encoding[END_REF][START_REF] Hofstetter | Magnetic resonance shear wave elastography using transient acoustic radiation force excitations and sinusoidal displacement encoding[END_REF] can provide an extended spatial mapping of elasticity. Another study, based on the generation of a "supersonic" series of ARFI pushes on a line, aimed at mapping tissue mechanical properties in a plane parallel to the ultrasound axis [START_REF] Liu | Supersonic transient magnetic resonance elastography for quantitative assessment of tissue elasticity[END_REF]. Both methods rely on the resolution of a complex inverse problem and the use of algorithms to estimate the shear wave velocity with several wave-propagation time steps.

In comparison, our method has been adapted for local estimation of viscoelasticity only within the HIFU ablation area i.e. directly where the ARFI push is applied. Although it does not allow for spatial mapping of the mechanical properties, our approach combines several advantages, such as the requirement of only one wave-propagation time step during the monitoring, which renders it particularly interesting in terms of temporal resolution, and the possibility to take tissue viscoelasticity into account through the evaluation of the τ parameter. Under the assumption of a linear viscoelastic model, this time constant is considered as an intrinsic property of the tissue, equal to the ratio between dynamic viscosity η and shear modulus μ. As a consequence, the value of τ is the same for both loading (ultrasound push ON) and relaxation (ultrasound OFF) phases, as already modelled in prior work [START_REF] Souchon | Transient MR elastography (t-MRE) using ultrasound radiation force: Theory, safety, and initial experiments in vitro[END_REF].

A major limitation of our model is that it does not consider shear wave dispersion, despite the fact that tissue is assumed to exhibit a viscoelastic response. The transient ARFI excitation leads to the propagation of shear waves with a given frequency content that was not considered in our model. Indeed, the shear wave propagation was modelled as a simple temporal shift of the focal displacement (Eq. 3-4). The major effect of this approximation is that our model ignores the distortion of the temporal displacement profile as it travels away from the focus, since it is assumed to always follow a strict first order exponential behavior (Eq. 3-4). An important

Fig. 1

 1 Fig. 1 General concepts of MR-ARFI

Fig. 2

 2 Fig. 2 Identification process of the viscoelastic properties using the proposed MR-ARFI-based method

  5 T MRI systems (MAGNETOM Aeraand MAGNETOM Sola, Siemens Healthineers, Germany). Main acquisition parameters were: TR 500 ms, TE 28 ms, FOV 256 mm × 256 mm (320 mm x 320 mm for the in vivo experiment), matrix 128 × 128, slice thickness 5 mm (6 mm for the in vivo experiment), EPI factor 37, MSG: G 36 mT/m -TMSG 10 ms -slice encoding direction. The spine coil was used in all the experiments. To improve the SNR, it was combined alternatively with a 19-cm flexible loop coil for both phantom experiments, or with the body matrix coil for the in vivo experiment.At every time step, two MR-ARFI images are sequentially acquired with opposite MSG polarity. This feature has two advantages: (1) it allows doubling the displacement sensitivity in MR-ARFI, by performing the difference between the two phase images and (2) it allows retrieving simultaneously temperature maps through proton resonance frequency shift (PRFS)

Fig. 4

 4 Fig. 4 Calibrated phantom experiment (Top row) 2D MR-ARFI spots acquired during the 2 nd repetition for 4 increasing values of toffset

  0  0.2 8.4  0.3 8.7  0.3 8.5  0.4 7.8  0.3 7.5  0.2 8.0  0.3 9.0  0.4 Rep. 2 8.0  0.2 8.5  0.3 8.7  0.3 8.5  0.4 7.8  0.3 7.5  0.2 7.9  0.3 9.1  0.4 Rep. 3 8.0  0.3 8.5  0.3 8.8  0.4 8.5  0.4 7.9  0.3 7.5  0.3 7.8  0.4 8.9  0.5 Rep. 4 8.0  0.2 8.5  0.3 8.8  0.4 8.5  0.5 7.8  0.4 7.5  0,3 7.8  0.3 8.9  0.6 Global average (kPa) 8.0  0.2 8.5  0.3 8.8  0.3 8.5  0.4 7.8  0.3 7.5  0.3 7.9  0.3 9.0  0.5

Fig. 5

 5 Fig. 5 HIFU ablation in a gelatin phantom

Fig. 6

 6 Fig. 6 HIFU ablation in ex vivo porcine muscle

Fig. 7

 7 Fig. 7 In vivo HIFU ablation

Fig. 8

 8 Fig. 8 Temporal evolution of the temperature and shear modulus μ() at the focal spot

Table 1 .

 1 The toffset parameter was manually chosen for each experiment according to the extent of the MR-ARFI pattern for the given Tpush/TMSG configuration.

		Phantom experiment	Ex vivo experiment	In vivo experiment
	ARFI toffset	-2 ms	1 ms	-3 ms
	TR		500 ms
	ARFI shots		Tpush 5 ms (duty cycle 1%) 200 W acoustic power
	Pause	70 ms	50 ms	50 ms
	HIFU heating	Theating 370 ms (duty cycle 74%) 30 W acoustic power	Theating 300 ms (duty cycle 60%) 80 W acoustic power	Theating 300 ms (duty cycle 60%) 120 W acoustic power

FUS sequence for the monitoring experiments (note that a qualitative time scale is used here for readability reasons) 12/27 HIFU and ARFI sequence parameters are provided in

Table 1 .

 1 HIFU and ARFI sequence parametersPolar maps of the shear modulus μ() were drawn for 64 angles equally distributed in the

	interval   [0°: 360°[.

Table 2 .

 2 Mean ± standard deviation for each toffset of the shear modulus (kPa) found with MR-ARFI in the commercial phantom, assuming isotropy (average over all radial directions), for each of the 4 repetitions and global values over the 4 repetitions

	15/27

This work has benefitted from funding of the FUI (Fonds Unique Interministériel, BPI France) for the UFOGUIDE project, and the ANR (Agence Nationale de la Recherche) French national program "Investissements d'Avenir" for the LABEX-CAMI (ANR-11-LABX-0004), the IHU Strasbourg (Institute of Image Guided Surgery, ANR-10-IAHU-0002) and the LabCom TechnoFUS lab (ANR-21-LCV3-0007-01).

improvement would be to refine the shear wave propagation model, by using for instance simulations based on Green's functions [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] in order to take shear wave dispersion into account.

One of the major side findings of this study is that both the viscoelastic mechanical relaxation time τ and the dynamic viscosity η, defined as η = τ.μ, were found to vary during the ablation in both ex vivo and in vivo porcine muscle. Although this is a promising finding, since it suggests that viscoelastic parameters may also serve as biomarkers of tissue integrity, it represents a significant limit to our method, which assumes that  is steady during the ablation.

Identifying both  and μ at the same time during the ablation seems therefore a necessary improvement to the proposed method. However, it may lead to much longer acquisitions, thus not compatible with the real-time requirements of HIFU ablation monitoring. Indeed, the current MR-ARFI-based method for estimating  relies on the use of numerous acquisitions with different toffset values [START_REF] Kaye | Adapting MRI Acoustic Radiation Force Imaging For In Vivo Human Brain Focused Ultrasound Applications[END_REF][START_REF] Dadakova | Optimization of acoustic radiation force imaging: Influence of timing parameters on sensitivity[END_REF]. From an optimization point of view, fitting these 2 parameters remains challenging on a single 1D-MR-ARFI phase profile, given the limited spatial resolution and experimental phase-to-noise ratio of the MR images. For this purpose, reshaping the identification process should be considered, which is beyond the scope of this paper.

In conclusion, a fast, non-invasive method has been proposed for the monitoring of tissue shear modulus and temperature during HIFU ablations. It provides quantitative shear modulus polar maps with every pair of MR-ARFI phase images acquired with opposite MR-ARFI-encoding. This fast local elasticity characterisation approach allows for monitoring tissue anisotropy, as illustrated in a tissue-mimicking gel phantom and in ex vivo and in vivo muscle porcine model.

Compliance with ethical standards

The in vivo experiment was led following the 3R principle for more humane research, and in following local legal ethics regulation (registered project APAFIS #14092-2018031513247711 v1). The animal was maintained under gaseous anaesthesia during imaging (1.8% isoflurane), with appropriate pain management (Propofol 1mg/kg) during HIFU. At the end of the experiment, the animal was euthanized through anaesthetics overdose (10 min, 5% isoflurane) coupled to IV injection of saturated KCl solution (0.5 mL/kg).