Cesare Corrado

Jean-Frédéric Gerbeau

Philippe Moireau

Identification of weakly coupled multiphysics problems. Application to the inverse problem of electrocardiography

Identification of weakly coupled multiphysics problems.

Application to the inverse problem of electrocardiography

Introduction

In the last few years, more and more attention has been paid to the problem of state and parameters identification for complex three-dimensional models used in biomedical applications. Several works can be cited: for example in cardiac electrophysiology [START_REF] Boulakia | Reduced-order modeling for cardiac electrophysiology. application to parameter identification[END_REF][START_REF] Burger | Stability analysis of the inverse transmembrane potential problem in electrocardiography[END_REF][START_REF] Chapelle | A Galerkin strategy with Proper Orthogonal Decomposition for parameter-dependent problems -Analysis, assessments and applications to parameter estimation[END_REF][START_REF] Nielsen | On the use of the bidomain equations for computing the transmembrane potential throughout the heart wall: An inverse problem[END_REF] or in cardiac mechanics [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-mri using a biomechanical heart model[END_REF][START_REF] Krishnamurthy | Patient-specific models of cardiac biomechanics[END_REF][START_REF] Xi | Myocardial transversely isotropic material parameter estimation from in-silico measurements based on a reduced-order unscented Kalman filter[END_REF], or also in hemodynamics [START_REF] Bertoglio | Sequential parameter estimation for fluid-structure problems. application to hemodynamics[END_REF][START_REF] D'elia | Applications of variational data assimilation in computational hemodynamics[END_REF][START_REF] Moireau | Sequential identification of boundary support parameters in a fluid-structure vascular model using patient image data[END_REF][START_REF] Perego | A variational approach for estimating the compliance of the cardiovascular tissue: An inverse fluid-structure interaction problem[END_REF]. The observations available in these contexts are often multiphysics since several modalities can be used simultaneously: electrocardiograms, electrograms, MRI, CT scan, flow measurements with ultrasound, pressure measurement with catheters, myocardium thickness measurement with piezoelectric sensors, etc. Up to now, the multiphysics nature of these problems has sometimes been taken into account in the direct models but it has rarely been used in the inverse problems.

The first purpose of the present study is to exhibit an example where multiphysics observations actually improve the identifiability of a coupled system. More precisely, considering an electromechanical model of the heart, it is shown that the estimation of an electrical parameter is improved if the electrical observation, namely the electrocardiogram, is enriched with mechanical observations, namely the movement of the myocardium. In this model, the electrophysiology acts as an input for the mechanics, but the electromechanical feedback is neglected [START_REF] Chapelle | An energy-preserving muscle tissue model: formulation and compatible discretizations[END_REF]. The coupling is therefore only one-way. The second purpose of this article is to show that a reduced order filtering strategy is well-suited to this class of multiphysics problems. Optimal filtering, like Kalman filter and its nonlinear extensions, is known to be very efficient but also too expensive to be used for large problems like those considered here. An effective strategy consists in using this kind of methods for the parameters only, and to address the uncertainties on the state variables through a less expensive approach. This strategy has been successfully used for example in mechanics where the state variable was handled with a Luenberger filter [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF]. In this paper, a similar strategy is proposed to filter the electrical state variables, by exploiting the special structure of the one-way coupled problem. But due to the structure of the equations, a Luenberger approach is not straightforward in electrophysiology. A different reduced order method, based on Proper Orthogonal Decomposition is then proposed, as it was done for cardiac mechanics in [START_REF] Chapelle | A Galerkin strategy with Proper Orthogonal Decomposition for parameter-dependent problems -Analysis, assessments and applications to parameter estimation[END_REF].

Background and related work

The goal of the inverse problem of electrocardiography, also called cardiac electrical imaging, is to reconstruct the electrical activity of the heart from body surface potential maps. Various strategies have been proposed since four decades. All of them assume that measurements of the electrical potential u T are available on parts of the torso boundary ∂Ω T .

The different strategies can be distinguished by the cardiac electrical source models they rely on. One of the first approaches was to estimate equivalent electrical dipoles [START_REF] Gulrajani | Moving dipole inverse ecg and eeg solutions[END_REF][START_REF] Mirvis | Detection and localization of multiple epicardial electrical generators by a two-dipole ranging technique[END_REF]. Another popular approach is to estimate the heart surface potential, usually called epicardial potential (even though pericardial potential would be more appropriate as noted in [START_REF] Messnarz | A comparison of noninvasive reconstruction of epicardial versus transmembrane potentials in consideration of the null space[END_REF]). The potential u T within the torso Ω T is assumed to be solution of the Poisson problem:

∇ • (D T • ∇ u T ) = 0, in Ω T u T = u e , on ∂Ω H 0 (1)
where ∂Ω H 0 denotes the boundary of the heart and D T is the electrical conductivity of the torso. The inverse problem then consists in estimating u e on ∂Ω H 0 (see e.g. [START_REF] Barr | Relating epicardial to body surface potential distributions by means of transfer coefficients based on geometry measurements[END_REF][START_REF] Burnes | A noninvasive imaging modality for cardiac arrhythmias[END_REF][START_REF] Wang | Finite-element-based discretization and regularization strategies for 3-d inverse electrocardiography[END_REF]). This problem being notoriously ill-posed, various regularizations have been proposed: Tikhonov [START_REF] Messinger-Rapport | Regularization of the inverse problem of electrocardiology: A model study[END_REF], the use of temporal information [START_REF] Greensite | An improved method for estimating epicardial potentials from the body surface[END_REF][START_REF] Oster | The use of temporal information in the regularization of the inverse problem of electrocardiography[END_REF], truncated Singular Value Decomposition or truncated Total Least Square [START_REF] Pullan | Mathematically modelling the electrical activity of the heart. From cell to body surface and back again[END_REF].

Another approach to address the inverse problem of electrocardiography is to consider the following equation within the heart Ω H 0 :

-∇ • ((D i + D e ) • ∇ u e ) = ∇ • (D i • ∇ v m ), in Ω H 0 ,
which is one of two equations of the bidomain model (see Section 2.1). Then, instead of estimating the epicardial potential u e on the surface ∂Ω H 0 , the goal is to estimate the transmembrane potential v m within the heart. In [START_REF] Messnarz | A comparison of noninvasive reconstruction of epicardial versus transmembrane potentials in consideration of the null space[END_REF], the two approaches were compared: investigating the null-space of the inverse problem, the authors concluded that the transmembrane potential-based formulation is more promising because it is based on a stronger biophysical a priori. As the epicardial potential approach, the transmembrane potential approach is ill-posed and must be regularized. In [START_REF] Liu | Noninvasive reconstruction of three-dimensional ventricular activation sequence from the inverse solution of distributed equivalent current density[END_REF], four variants of L 2 -Tikhonov regularization are compared. In [START_REF] Nielsen | On the possibility for computing the transmembrane potential in the heart with a one shot method: An inverse problem[END_REF], H 1 -Tikhonov regularization is used with a prior taking two different homogeneous values in the myocardium depending on cardiac phase (plateau or rest values). The inverse problem, formulated in a PDE-constrained framework, is addressed by directly solving the optimality saddle-point problem. In [START_REF] Nielsen | On the use of the resting potential and level set methods for identifying ischemic heart disease: An inverse problem[END_REF], the estimation of the transmembrane potential is combined with a level set technique to efficiently identify the location of a myocardial infarction. In [START_REF] Wang | Inverse electrocardiographic source localization of ischemia: An optimization framework and finite element solution[END_REF], the approach of [START_REF] Nielsen | On the possibility for computing the transmembrane potential in the heart with a one shot method: An inverse problem[END_REF] is generalized to more general objective functions and constraints in order to identify ischemic regions (characterized by lower amplitude during the plateau phase). Two different regularizations are investigated: the Tikhonov regularization, which is found to overestimate ischemic regions but with good sensitivity, and the Total Variation regularization which is found to underestimate ischemic regions but with high specificity.

In [START_REF] Rincon | Computing the electrical activity of the heart with a dynamic inverse monodomain operator[END_REF], the authors note that the usual regularization techniques have no physical ground. Instead, they propose to regularize the inverse solution with the monodomain equations coupled to the Fenton-Karma ionic model. The strategy proposed in the present paper has some similarities with this approach. We rely on a full electrophysiological model of the action potential coupled to the Poisson problem [START_REF] Auroux | Back and forth nudging algorithm for data assimilation problems[END_REF] to estimate the solution of the heart electrical activation. This physical model is personalized on the fly with respect to its parameters in order to adapt it to a specific patient. Furthermore, we propose an additional step of modeling by considering the mechanical response to the electrophysiological activation, so we are able to also integrate mechanical measurements. Indeed, we believe that multimodal observations improve the identifiability of the complete model and therefore improve the quality of the electrical and mechanical state reconstruction. Another originality of our work is the use of a sequential data assimilation strategy that is adapted to a coupled electromechanical evolution model. Here we demonstrate how state-of-the-art gain filter on the electrophysiological model and on the mechanical model can be aggregated to propose a joint gain filter for the coupled problem.

Organization of the present work

The paper is organized as follows. In Section 2, the electromechanical model is presented. In Section 3, the observation operators -namely the measurements -are detailed for the electrical and the mechanical variables. In Section 4, the general notions of data assimilation, optimal filtering and reduced order filtering are reviewed. Although the algorithms of this section are not new, their presentation differs from what is most often done in the literature since a purely deterministic description is adopted. In Section 5, the algorithms used for the electromechanical problem are proposed and analyzed. In Section 6, numerical experiments based on synthetic data are presented. The main purpose is to estimate a non-homogeneous parameter of the electrical model using electrical and mechanical observations.

Models

We present the models in a time and space continuous context before entering into the discretization details and their numerical implementation. Concerning the continuous context, we denote by an underline character any vector field of R 3 and two underlines any second-order tensor.

Electrophysiology

A widely accepted model of the macroscopic electrical activity of the heart is the socalled bidomain model [START_REF] Colli Franzone | Degenerate evolution systems modeling the cardiac electric field at micro-and macroscopic level. evolution equations, semigroups and functional analysis[END_REF][START_REF] Sachse | Computational Cardiology: Modeling of Anatomy, Electrophysiology, and Mechanics[END_REF][START_REF] Sundnes | Computing the electrical activity in the heart[END_REF][START_REF] Tung | A bi-domain model for describing ischemic myocardial D-C potentials[END_REF]. It consists of two degenerate parabolic reactiondiffusion PDEs which describe the dynamics of the averaged intra-and extracellular potentials u i and u e , coupled to a system of ODEs defining an ionic model. This model is related to the chemical dynamics of the myocardium cell membrane, in terms of the (vector or scalar) variable w representing the distributed ion concentrations and gating states, or a phenomenological counterpart. The model reads

       a m c m ∂ t v m + I ion (v m , w) -∇ • (D i ∇ u i ) = a m I app , in Ω H 0 , a m c m ∂ t v m + I ion (v m , w) + ∇ • (D e ∇ u e ) = a m I app , in Ω H 0 , ∂ t w + g(v m , w) = 0, in Ω H 0 , (2)
where v m = u i -u e represents the transmembrane potential, c m is the membrane capacitance per unit area, a m is a constant representing the rate of membrane area per unit volume, D i , D e are the intra-and extra-cellular conductivity tensors, I app is an external volume current. In (2) the function g represents an ionic model. In this article, the Mitchell-Schaeffer ionic model [START_REF] Mitchell | A two-current model for the dynamics of cardiac membrane[END_REF] is considered, with the same rescaling as in [START_REF] Boulakia | Mathematical Modeling of Electrocardiograms: A Numerical Study[END_REF]. It is a reduced complexity model capable of integrating relevant phenomenological properties of the ventricle cell membrane:

               I ion (v m , w) = - w τ in (v m -V min ) 2 (V max -v m ) V max -V min + 1 τ out v m -V min V max -V min , g(v m , w) =      w τ open - 1 τ open (V max -V min ) 2 if v m ≤ V gate , w τ close if v m > V gate , (3)
where V gate , τ in , τ out , τ open , τ close are given constants and V min and V max are scaling constants (typically V min = -80mV and V max = 20mV ).

On the boundary ∂Ω H 0 , we have D i • ∇ u i • n = 0, and the heart is assumed to be isolated, D e • ∇ u e • n = 0, as often done in the literature [START_REF] Clements | Activation dynamics in anisotropic cardiac tissue via decoupling[END_REF][START_REF] Potse | A comparison of monodomain and bidomain reaction-diffusion models for action potential propagation in the human heart[END_REF]. For well-posedness, the condition Ω H 0 u e = 0 is enforced.

Summing and subtracting the first two equations of (2), the system reads

                     a m c m ∂ t v m + I ion (v m , w) -∇ • (D i • ∇ (u e + v m )) = a m I app , in Ω H 0 , -∇ • ((D i + D e ) • ∇ u e ) -∇ • (D i • ∇ v m ) = 0, in Ω H 0 , ∂ t w + g(v m , w) = 0, in Ω H 0 , (D e • ∇ u e ) • n = 0, on ∂Ω H 0 , (D i • ∇ (u e + v m )) • n = 0, on ∂Ω H 0 . (4)
A P 1 finite element discretization of the potential variables is used -see the corresponding refined electrical mesh in Figure 2. This leads to a discretization space V e h with N e degrees of freedom (i.e. V e h R N e ). To each field -for instance v m -is associated its corresponding approximation, denoted with an index h -for example v mh . Equivalently, this approximation can be represented by its corresponding vector of degrees of freedom written with a vector in uppercase letteri.e. V m . The finite element vector of degrees of freedom or the discretization of linear form are then defined with vectors in straight uppercase letter whereas the finite element matrix operator with a bold uppercase letter. With this notation:

∀ U ∈ V e h , U M e V m = Ω0 v mh u h dΩ, U M e a V m = Ω0 a m c m v mh u h dΩ, ∀ U ∈ V e h , U K e i V m = Ω0 D i • ∇ v mh • ∇ u h dΩ, ∀ U ∈ V e h , Ω0 a m (I app -I ion (v mh , w h )) u h dΩ U M e I I app -I ion ( V m , W )
and if a spatial discretization of the internal variable is done by node

∀ U ∈ V e h , Ω0 g(v mh , w h )u dΩ U M e G( V m , W )
where the applied current and the ionic variables have been interpolated, allowing to define the ionic variables at the nodes instead of the quadrature points. Finally, after spatial discretization system (4) reads

       M e a ˙ V m + K e i ( V m + U e ) = M e I ( I app -I ion ( V m , W )), (K e i + K e e ) U e + K e i V m = 0, ˙ W + G( V m , W ) = 0. (5)
The state of this system is X e = Vm W while U e appears as an auxiliary variable verifying the static equilibrium (5) 2 with V m . This can be summarized by

     Ẋe = A e (X e , θ e ), X e (0) = X e + ζ X e , θ e = θ e + ζ θ e . ( 6
)
where θ e denotes the vector of parameters of the electrophysiology model, ζ θ e and ζ X e the uncertainty on the parameters and the initial condition respectively.

Mechanics

The heart domain is denoted by Ω H (t) at any time t. This domain is the image of a reference configuration Ω H 0 through the solid deformation mapping φ φ :

Ω H 0 × [0, T ] -→ Ω H (t), (ξ, t) -→ x = ξ + y(ξ, t)
where y is the solid displacement. The solid velocity is given by v = ẏ. The deformation gradient F is given by F (ξ, t) = ∇ ξ φ = 1 + ∇ ξ y, and its determinant is denoted by J. The right Cauchy-Green deformation tensor is defined by C = F • F , the Green-Lagrange tensor by e = 1 2 (C -1) = 1 2 ∇ ξ y + (∇ ξ y) + (∇ ξ y) • ∇ ξ y , and its linearization by ε = 1 2 ∇ ξ y + (∇ ξ y) . The mass per unit volume is denoted by ρ and the Cauchy stress tensor by σ. In the reference configuration, the first and second Piola stress tensor are respectively defined by

T = Jσ • F -and Σ = F -1 • T = JF -1 • σ • F -.
The constitutive law is assumed to be a combination of a hyperelastic law of potential W , a viscous component chosen proportional to the strain rate ė, and an active part along the fiber direction τ represented by 3 internal variables which are e c the active strain, k c the active stiffness and τ c the associated active stress [START_REF] Chapelle | An energy-preserving muscle tissue model: formulation and compatible discretizations[END_REF]:

Σ(e, e c , k c , τ c ) = ∂W ∂e (e) + η s ė + σ 1D (e c , k c , τ c )τ ⊗ τ , (7)
with σ 1D = 1+2ec 1+2τ •e•τ (τ c + µ ėc ), where these 3 internal variables rely on a chemicallycontrolled constitutive law describing the myofibre mechanics [START_REF] Bestel | A biomechanical model of muscle contraction[END_REF][START_REF] Chapelle | An energy-preserving muscle tissue model: formulation and compatible discretizations[END_REF][START_REF] Huxley | Muscle structure and theories of contraction[END_REF]:

∂ t k c = -(|u| + α | ėc |)k c + n 0 k 0 |u| + in Ω H 0 ∂ t τ c = -(|u| + α | ėc |)τ c + ėc k c + n 0 σ 0 |u| + in Ω H 0 (8)
with α, k 0 , σ 0 given parameters, n 0 a function of e c accounting for the Frank-Starling effect and u directly related to the electrical activity of the heart by

u(t) = av m (t) + b
where a and b are two scaling parameters.

Concerning the boundary conditions, the external organs are modeled by visco-elastic boundary conditions on a sub-part of the epicardium: T • n = k s y + c s v on Γ n (t). A uniform pressure is enforced on the left and right endocardium: σ •n t = p v,i n t on Γ n,i (t), i = {1, 2}. In summary, the mechanical problem reads

                       ∂ t y = v, in Ω H 0 ρ∂ t v -∇ • (T ) = 0, in Ω H 0 ∂ t k c = -(|u| + α | ėc |)k c + n 0 k 0 |u| + , in Ω H 0 ∂ t τ c = -(|u| + α | ėc |)τ c + ėc k c + n 0 σ 0 |u| + , in Ω H 0 T • n = k s y + c s v, on Γ n T • n = Jp v,i F -• n, on Γ c,i T • n = 0, on ∂Ω H 0 \((∪ i Γ c,i ) ∪ Γ n ) (9)
with the constitutive law [START_REF] Bestel | A biomechanical model of muscle contraction[END_REF].

Remark 1

Several improvements can be formulated on this model: other active or passive constitutive laws, more sophisticated boundary conditions for the tethering of the myocardium or hemodynamics, etc. The estimation strategy presented here can be adapted to any of these improvements since we only refer to state-space model description to formulate our data assimilation methods.

The system is discretized with a P 1 finite element -see the corresponding mechanical mesh in Figure 2 -with a 5% compressibility acceptance in order to avoid any numerical locking. Concerning the fibre directions, we prescribed them on each point of the mesh with an elevation angle varying from -60 degrees to 60 degrees through the myocardium thickness. The discrete system is based on the variational formulation associated with [START_REF] Bonnabel | Observer-based Hamiltonian identification for quantum systems[END_REF]:

∀v ∈ V v , Ω0 ρ∂ t v • v dΩ + Ω H 0 Σ(e, ė, e c , k c , τ c ) : d y e • v dΩ + Γn (k s y + c s v) • v dΓ = - i Γc,i Jp v (F -• n) • v dΓ (10)
Using the same convention as for the electrophysiological model discretization, the mass operator is defined by

∀ V , V M m ˙ V = Ω H 0 ρ∂ t v h • v h dΩ,
the stress residual by

∀ V , V K m ( Y , V ) = Ω H 0 Σ(e h , ėh , e c,h , k c,h , τ c,h ) : d y e h • v h dΩ,
some weighted mass operators on the boundary by

∀ V , V M m ks,Γn Y = Γn k s y h • v h dΓ, V M m cs,Γn V = Γn c s v h • v h dΓ,
and a following pressure operator by

∀ V , V N m ( Y ) = Γc,i J h (F - h • n) • v h dΓ.
The internal variables e c , k c , τ c are gathered in a vector ı c , which is discretized at the integration points. A vector of degrees of freedom ı c is associated with the discontinuous field (e c,h , k c,h , τ c,h ) and the model ( 8) is discretized by ˙

ı c = B m ( ı c , Y , V ). The complete mechanical model is spatially discretized into        ˙ Y = V , M m ˙ V + K m ( Y , V , ı c ) + M m ks,Γn Y + M m cs,Γn V = -i p v,i N m ( Y ), ˙ ı c = B m ( ı c , Y , V ). ( 11
)
The state of this system is

X m = Y V ıc
. Denoting by θ m the set of the parameters characterising the mechanics, affected by an a priori uncertainty ζ θ m and with ζ X m the a priori uncertainty on the initial condition, the mechanical system reads

     Ẋm = A m (X m , X e , θ m ), X m (0) = X m + ζ X m , θ m = θ m + ζ θ m , (12)
where the electrical variable X e , solution to [START_REF] Bertoglio | State observers of a vascular fluid-structure interaction model through measurements in the solid[END_REF], can be seen as an input.

Electromechanical coupling

From a computational point of view, the electromechanical problem consists in solving two coupled systems. We chose to keep the two sub-systems in independent solvers. This choice allows us to use legacy codes, to make their maintenance easier and to take advantage of the specific numerical methods adapted to each physical compartment. The coupling algorithm sketched in Figure 1 is handled by a "master" code which exchanges the heart displacements and the transmembrane potential with the electrical and mechanical software. In this work, it is assumed that there is no electromechanical feedback. The transmembrane potential v m is sent to the mechanical problem. Then, the one-way coupling is performed through the quantity u = av m + b which triggers the mechanical contraction via a change in the active stiffness k c and in the active stress τ c , see [START_REF] Blum | Data assimilation for geophysical fluids[END_REF] Because of the characteristics of the action potential propagation, the space and time steps required by the electrophysiology are typically smaller than those needed by the mechanics. To ensure accuracy at a reasonable computational cost, each sub-problem is solved with its own space and time step. The data are transmitted by the master code at some check-point in time (see Figure 1) and interpolated in space from a mesh to another (see Figure 2). Eventually, a typical complete direct simulation of our model is presented in Figure 3. The inverse problem will be addressed on the discrete formulations. The state and the parameters are defined as the combination of the electrical and the mechanical ones

X = X e X m , θ = θ e θ m . Defining A(X, θ) = A e (X e ,θ e )
A m (X m ,X e ,θ m ) , this state variable follows the dynamics

Ẋ = A(X, θ) (13)
The initial condition includes the quantities ζ X and ζ θ which model the uncertainties on the initial condition and on the parameters respectively. When necessary, a trajectory will be denoted by

X [ζ X ,ζ θ ] .
This system can also be written in an augmented form by concatenated the state and parameters in X = X θ . This augmented state satisfies Ẋ = A(X). The real trajectory is assumed to be a solution of the model for a given value of the uncertainties and is given by X † = X † θ † .

Measurements

Electrophysiological measurements

The best way to access the electrophysiological quantities would be to directly measure the electrophysiological potential at the heart surface. This is possible with an invasive procedure where a basket of electrodes is immersed in the ventricle and records the potential close to the endocardium surface ([21] e.g.). In that case, the electrophysiological measurements can be represented by

Z e (t) = H e x e + χ e (14)
where the observation operator H e interpolates the values at the basket electrodes and χ e gathers all the measurement errors. After spatial discretization

Z e (t) = H e X e + χ e (15)
where now χ e takes also into account the finite element discretization error.

To avoid such invasive measurement procedures, it is preferable to rely on electrocardiograms. Physiologically, the electrical potential diffuses from the heart to the rest of the body through the pericardium. The ECG depicts the time course of standardized potential differences on the body surface. We refer to [START_REF] Boulakia | Mathematical Modeling of Electrocardiograms: A Numerical Study[END_REF] for examples of healthy and pathological ECG obtained with the electrophysiology model used in the present paper.

The electrical diffusion within the body Ω T can be modelled by a Poisson equation with the extracellular potential as a boundary condition on the heart surface and homogeneous Neumann boundary elsewhere:

     ∇ • (D T • ∇ u T ) = 0, in Ω T u T = u e , on ∂Ω H 0 (D T • ∇ u T ) • n = 0, on ∂Ω T \∂Ω H 0 (16)
where the conductivity tensor is assumed to be isotropic, but non-homogeneous to account for the different conductivities of the lungs, the bones and the rest of the body (see [START_REF] Boulakia | Mathematical Modeling of Electrocardiograms: A Numerical Study[END_REF]).

After spatial discretization -see the corresponding thorax mesh in Figure 2 -the ECG measurements denoted by Z e are still related to the electrical state variable by an observation operator H e which is the composition of the discrete diffusion operator and a linear combination of the potential at 9 points of the body surface. Assuming that the model is accurate enough, the actual ECG can still be represented by [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-mri using a biomechanical heart model[END_REF] where χ e gathers all the possible measurement errors, including the errors resulting from the modelling of the measurement procedure. After spatial discretization, a relation of the form (15) still holds, but with χ e also accounting for the discretization error.

Mechanical measurements

The heart contraction can also be perceived from a kinematic and even a mechanical perspective. In fact, the displacements of the cavity can be observed using imaging modalities, for example Cine-MRI sequence, CT-sequence or even Tagged-MRI sequence. Considering for example this last type of measurement and assuming adequate registration strategies allowing to extract 3D displacements of material points [START_REF] Rutz | Accelerated whole-heart 3D CSPAMM for myocardial motion quantification[END_REF], we can assume that displacements are measured in a large part of the ventricles, typically between two short axis planes, the first one slightly lower than the base and the second one slightly higher than the apex -see Figure 3. With more standard sequences, like Cine-MRI or CT sequences, we can benefit from optical flow strategies to extract displacements most of the time on part of the boundary. In both cases, after discretization, we can assume the existence of an observation operator H m such that the measured displacements are given by

Z m (t) = H m X m + χ m ,
where χ m is now a vector taking into account the discretization error.

Remark 2 (Generalized observation operators)

Extracting displacements using registration techniques is still a challenge and it often leads to the measurement of an apparent motion instead of the real motion. As a consequence, it can be more adequate to consider the structure measured in itself more than its collection of material points. This was done typically in [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF] where in Cine-MRI and CT sequences the contours of interest are extracted and compared to the contours produced by the simulation by the use of a discrepancy field of the form

D(ξ, t) = dist St (x(ξ, t), t) = χ(ξ, t).
Hence, the closer to 0 the discrepancy is the more information we have extracted from the images. As demonstrated in [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF], this type of discrepancy operator is a simple non-linear generalization of the discrepancy computed by subtracting to the actual measurements the one produced by the use of the observation operator defined below. This generalization has demonstrated its efficiency in several real cases of data assimilation investigations [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-mri using a biomechanical heart model[END_REF][START_REF] Moireau | Sequential identification of boundary support parameters in a fluid-structure vascular model using patient image data[END_REF]. But, for the sake of simplicity, and since this work is primarily aimed at providing a proof of concept, we limit ourselves to more simple observation operators. Note that the use of a non-linear observation operator is, in the end, always analyzed in the light of the corresponding linearized operator [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF].

The stresses are more intricate to measure. For now, the most common measurement available is an intracavity pressure obtained after invasive catheterization or reconstructed from arterial blood pressure measurements. In this article we assume that we can consider the intracavity pressure as a given source term in our model [START_REF] Bonnabel | Observer-based Hamiltonian identification for quantum systems[END_REF]. In a more general case, we could model the evolution of the intracavity pressures with the help for example of a Windkessel model and then use the possibly noisy measurements in a data assimilation context where the intracavity pressures are additional variables of the model.

Multi-modalities measurements

Finally the observations are concatenated and therefore given and the complete form

Z = Z e Z m = H e 0 0 H m X e † X m † + χ e χ m = HX † + χ, (17)
which for the sake of generality will not be necessarily considered as linear when it is not mandatory.

Data assimilation principles

Data assimilation has become a very popular strategy to estimate a wide range of modeling uncertainties in numerical simulations [START_REF] Blum | Data assimilation for geophysical fluids[END_REF]. It has been initiated for environmental sciences but has now reached life sciences and especially cardiology. To introduce its main concepts in the state-space formalism, we consider a dynamical model ẋ = A(x, ϑ, t).

(

) 18
In this equation, x denotes the state variable, namely, the physical quantity which the model aims at describing during its time evolution. In this generic notation, the whole model is essentially summarized in the dynamical operator A, which applies on the state variable itself, and may depend on time t as well as on a set of physical parameters denoted by ϑ. In this work, (18) corresponds to the electromechanical system defined by ( 4) and ( 9), or -to avoid issues associated with PDEs -its discrete counterpart ( 6) and ( 12) (but most of what is presented here can be actually generalized to infinite dimensional systems [START_REF] Bensoussan | Filtrage optimal des systèmes linéaires[END_REF]).

The initial condition x(0) and the parameter vector ϑ must be prescribed in [START_REF] Clements | Activation dynamics in anisotropic cardiac tissue via decoupling[END_REF]. Hence, when these quantities are unknown, they have to be estimated. In general x(0) and ϑ are decomposed into known parts x and ϑ called a priori and uncertain parts ζ x and ζ ϑ :

x(0) = x + ζ x , ϑ = ϑ + ζ ϑ .
When necessary, the trajectory of x will be denoted by

x [ζ x ,ζ ϑ ] .
Let us now re-consider an augmented state x = ( x ϑ ) in order to integrate parameter uncertainties. The corresponding augmented model is defined by

ẋ = A( x),
and

x(0) = x + ζ x ϑ + ζ ϑ = x + ζ. (19)
Hence the parameters and the state can be consider in the same formalism. However, these two quantities differ in their dimension. Indeed, when the state is defined based on a PDEs model, the dimension of the state variable after space discretization is typically of the size of 10 3 to 10 7 degrees of freedom. By contrast, the size of the parameter vector is generally much more limited. Even if distributed parameters are considered, their variation should be considered smooth enough so that they can be discretized on a coarse mesh or a subdivision of the domain into large regions. Therefore, typically less than a hundred of parameter values need to be estimated. The measurements presented in Section 3 can be cast in a general form: given a real trajectory x † , the noisy observations are represented using an observation operator H such that

z = H(x † , t) + χ(t), (20)
where z denotes the actual data field. Note that H can possibly be extended to a function of x and will then be denoted by H.

Data assimilation by filtering

A common strategy in data assimilation is to minimize a criterion that balances the confidence in the a priori value of the state and parameters and a discrepancy measure between the given observations and the simulated ones. In the state-space formalism this criterion is typically least-square and reads

J (ζ) = 1 2 ζ 2 P -1 + 1 2 T 0 z -H( x [ζ] ) 2 M dt. ( 21
)
where M is a metric on the observation space and P is the inverse of a metric which can be interpreted as an initial uncertainty covariance. The so-called 4D-Var method [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF] consists in finding the unknown quantities ζ by minimizing J under the constraint of following the dynamics [START_REF] Clements | Activation dynamics in anisotropic cardiac tissue via decoupling[END_REF]. When T goes to infinity the minimization is expected to produce a better and better estimate. This minimization under constraint requires the computation of an adjoint variable used to compute the gradient. Hence, a gradient based descent algorithm requires numerous iterations of the direct and adjoint dynamics. When we do not wish to precisely retrieve the initial error ζ but only seek to accurately approximate the state x = (x, ϑ) "independently" of the possible initial error, we can avoid minimization iterations by the use of sequential estimation methods. The principle in sequential estimation methods is to introduce a modified system -denoted by a circumflex accent -called observer whose dynamics is changed to incorporate a correction based on the measured discrepancy. The new system dynamics therefore reads

ẋ = A( x, t) + G(z -H( x)), x(0) = x ( 22
)
where G is called the observer gain or filter. The ultimate objective of the observer x(t) is to converge in time to the real trajectory

x(t) t→∞ ---→ x(t)
and the gain has to be designed regarding this objective. Two strategies are commonly followed in this respect. First, the gain can be constructed from the optimality criterion ( 21) by defining the so-called optimal observer1 -or optimal sequential estimator -with

x(t) = x [argmin J (•,t)] (t).
An optimal gain is obtained by differentiating this definition with respect to the time variable t appearing in both the trajectory x(t) and in the criterion J (•, t). This approach is well-known in a linear frameworki.e. where all the operators are linear -and leads to the famous Kalman-Bucy filter [START_REF] Bensoussan | Filtrage optimal des systèmes linéaires[END_REF][START_REF] Kalman | New results in linear filtering and prediction theory[END_REF][START_REF] Simon | Optimal State Estimation: Kalman, H ∞ , and Nonlinear Approaches[END_REF]. The optimal filter is then defined from the solution of a Riccati equation. In a non-linear framework, the gain is more intricate to compute and derives from a Hamilton-Jacobi-Bellman solution [START_REF] Fleming | Deterministic nonlinear filtering[END_REF]. Hence, in this case, numerous works rely on approximate solutions as defined for example by the Extended Kalman Filter (EKF) which uses the Riccati equation of the linear Kalman filter with the tangent operator of the non-linear model and observation operators [START_REF] Simon | Optimal State Estimation: Kalman, H ∞ , and Nonlinear Approaches[END_REF]. Eventually, the great advantage of the optimal sequential approach is that, like the 4D-Var, it can be defined for every model and every observation operator. The biggest drawback is that -even with approximate solutions -the filter is very costly to compute, especially for large dimensional systems like the ones produced by PDEs or their discretizations. Therefore, alternative strategies -the Luenberger observers and the Reduced Order Optimal filtering -will be presented below.

The Extended Kalman Filter (EKF)

Even if we will have to rely on alternative strategies, it is interesting to develop the optimal filtering equations. For the sake of simplicity, we proceed after space discretization. The space discretized system variable is denoted by X and the space discretized parameters by θ. The observer state X and parameters θ follow a modified version of the dynamics:

           Ẋ = A( X, θ) + G X (Z -H( X, t)), θ = G θ (Z -H( X, t)), X(0) = X , θ(0) = θ , (23)
where the gain G = G X G θ is a linear operator. In the augmented form with X = X θ we write

Ẋ = A( X) + G(Z -H( X, t)).
The most classical gain is given by the Extended Kalman Filter (EKF): G = P(dH) M, where P is obtained from the solution of a Riccati equation on the augmented form

Ṗ = (d A)P + P (d A) -P(d H) M(d H)P. ( 24
)
The operator P solution of the Riccati equation is called covariance since it can be linked in a stochastic framework to the covariance of the state estimation error evolving during the sequential estimation [START_REF] Simon | Optimal State Estimation: Kalman, H ∞ , and Nonlinear Approaches[END_REF]. When decomposed on the state and parameter

P = P XX P Xθ (P Xθ ) P θθ this gives G X = P XX (d X H) M, G θ = (P Xθ ) (d X H) M, and
         Ṗθθ = -(P Xθ ) (d X H) M(d X H)P Xθ , P θθ (0) = P * ṖXθ = (d X A)P Xθ + (d θ A)P θθ -P XX (d X H) M(d X H)P Xθ , P Xθ (0) = 0 ṖXX = (d X A)P XX + (d θ A)(P Xθ ) + P XX (d X A) + P Xθ (d θ A) -P XX (d X H) M(d X H)P Xθ , P XX (0) = P (25)
The practical algorithms -for instance time-discrete EKF or UKF for Unscented Kalman Filter -derived from this formulation are presented in Appendix A. Even if UKF differs from EKF, it is easy to prove that their analysis which relies on the linearization of their respective estimation error is based on the same error system [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF]. This is due to the fact that when all operators are assumed to be linear, a finite difference operator or the tangent operator are in fact identical. This implies that we will rely on UKF for the numerical simulation whereas, for the estimator analysis, we will keep EKF which can easily be written in both time-continuous and time-discrete formulations.

Various sequential estimators can be defined on any dynamical system, and parameter identification can be easily included. This is due to the underlying optimal principles behind the estimator presented before, even if in a non-linear framework we have only relied on approximation rules to extend the optimal filters found in the linear case. However, these estimators have a fundamental drawback since the covariance operator P is a full matrix of M N X +d (R). Recalling that N X is the dimension of the state variables which discretize the PDEs field variables, the optimal estimators are intractable for classical finite element models.

Reduced-Order Optimal Filters

A classical first strategy for circumventing the curse of dimensionality consists in assuming a specific reduced-order form for the covariance operators. For example, making the ansatz ∀t,

P(t) = L(t)U(t) -1 L(t) (26)
with U an invertible matrix of small size r and L an extension operator, we can show [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF] that with linear operators the solution of the Riccati equation ( 24) in augmented form reduces to L = AL and U = L H MHL, [START_REF] Huxley | Muscle structure and theories of contraction[END_REF] which are now computable in practice. In a non-linear framework, the covariance dynamics can then be approximated as in [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF] by extending [START_REF] Huxley | Muscle structure and theories of contraction[END_REF] as

L = (d A)L and U = L (dH) M(dH)L. ( 28
)
This observer is called Reduced-Order Extended Kalman Filter (RoEKF). The time discrete version of this approach is presented in Appendix B, as well as its UKF counterpart.

Luenberger Filters

Another way to circumvent the curse of dimensionality is to built a filter which is not based on an underlying optimal criterion. This idea was initially introduced in [START_REF] Luenberger | An introduction to observers[END_REF] and therefore is often called Luenberger filter and observer. It was quickly popularized in data assimilation [START_REF] Hoke | The initialization of numerical models by a dynamic-initialization technique(fluid flow models for wind forecasting)[END_REF][START_REF] Stauffer | Use of four-dimensional data assimilation in a limited-area mesoscale model. part i: Experiments with synoptic-scale data[END_REF][START_REF] Auroux | Back and forth nudging algorithm for data assimilation problems[END_REF] where the curse of dimensionality was very limiting for systems coming from the discretization of PDEs. In data assimilation this strategy is also referred to as the nudging approach because the filter is designed to "gently" correct the dynamics. The goal of the correction is to make the estimation error x = x † -x converge to 0. In a linear framework, the state error dynamics is ẋ = (A -GH) x, thus the strategy is to design a filter gain G based on the underlying physics of the problem which makes the dynamics operator (A -GH) dissipative -strongly if possible. Such Luenberger filters are designed for the state estimation of specific physical systems and for specific observation operators. For instance, filters have been proposed for transport equations [START_REF] Auroux | Back and forth nudging algorithm for data assimilation problems[END_REF], Schrödinger problems [START_REF] Bonnabel | Observer-based Hamiltonian identification for quantum systems[END_REF], for waves [START_REF] Chapelle | Exponential convergence of an observer based on partial field measurements for the wave equation[END_REF], beams, elasticity [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF][START_REF] Xu | Infinite dimensional observers for vibrating systems[END_REF] or fluid-structure problems [START_REF] Bertoglio | State observers of a vascular fluid-structure interaction model through measurements in the solid[END_REF]. In Section 5.2, this approach will be used for the mechanical problem.

4.5. Joint state and parameter estimation and coupled system estimation with Luenberger filters System [START_REF] Colli Franzone | Degenerate evolution systems modeling the cardiac electric field at micro-and macroscopic level. evolution equations, semigroups and functional analysis[END_REF] has been defined from the augmented state vector gathering the original state and the parameters. After decomposition on x and ϑ, this implies that the initial parameters now evolve in time within the observer

ẋ = A(x, θ, t) + G x (z -H(x)), x = x , θ = G ϑ (z -H(x)), θ(0) = ϑ , (29)
We have seen that the optimal filter does not make any difference between state and parameters. However, when dealing with Luenberger observer, there is no direct way to incorporate the parameter identification -i.e. to design G ϑ -from an already defined Luengerger filter G x on the state. This problem is known as adaptative filtering [START_REF] Baumeister | Adaptive methods for parameter identification[END_REF][START_REF] Zhang | Adaptive observer for multiple-input-multiple-output (MIMO) linear time-varying systems[END_REF].

In particular in [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF][START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF] a systematic strategy has been proposed to extend a possible Luengerger observer to a joint state and parameters estimation by combining the Luenberger filter and an optimal filter reduced to the remaining parameter space. Taking into account that the parameter space is of much smaller dimension, this strategy allows us to compute a physically-based gain on the large dimensional state and an optimal-based filter on the small parameter space. The resulting observer is proved to converge under adequate assumptions in linear cases and can be extended to non-linear cases.

When considering multi-physics coupled systems, as the electromechanical problem (13), the same type of problem occurs as for joint state and parameters estimation. If optimal filters are considered, there is no difficulty to combine them. When measurements are available, and when Luenberger filters have also been designed for the different types of physics, it is also easy to combined them. But, as said before, optimal filters are extremely expensive, and it is not always possible to design a Luenberger filter for all the physical compartments of the problem. It will be justified in Section 5.3 that, as done for the state and parameter estimation, it is possible for one-way coupled problems to combine a Luenberger filter in one part of the system and a reduced optimal filter on the other part.

Data assimilation for the electromechanical model

State estimation in electrophysiology: a reduced-order approach

The above data assimilation principles can been applied to the system modeling the cardiac electrophysiology. Using a few measures of the electrical potential, the goal is to reduce the uncertainties on the state and the parameters involved in the modeling of cardiac cells. For this system, contrary to the mechanical system that will be presented in the next section, there is no straightforward Luenberger observer. The state has therefore to be filtered differently.

A possible strategy to overcome the curse of dimensionality induced by the optimal filtering is to discretize the problem on a low-dimensional basis. Here, we propose to build this basis by Proper Orthogonal Decomposition (POD). The POD basis is obtained by keeping the most relevant modes resulting from a Principal Component Analysis of a set of pre-computed solutions. For more details about POD, we refer for example to [START_REF] Kunisch | Galerkin proper orthogonal decomposition methods for parabolic problems[END_REF][START_REF] Rathinam | A new look at proper orthogonal decomposition[END_REF] and to [START_REF] Boulakia | Reduced-order modeling for cardiac electrophysiology. application to parameter identification[END_REF][START_REF] Chapelle | A Galerkin strategy with Proper Orthogonal Decomposition for parameter-dependent problems -Analysis, assessments and applications to parameter estimation[END_REF] for applications to electrophysiology.

We denote by Π ∈ M N e ,NPOD (R) the matrix made of the first N POD modes. These modes are orthonormal with respect to a given scalar product, typically in l 2 (R NPOD ) or L 2 (Ω H ). Denoting by M POD the Gramian matrix associated with this scalar product, the POD expansion coefficients of a vector X e are given by α = Π M POD X e . Defining X e POD = Πα and X e ⊥ = X e -Πα, the state vector can be decomposed as X e = X e POD +X e ⊥ . As for the reduced filtering, a first strategy would consist in replacing the full order model by the reduced one:

     Ẋe POD = A e POD (X e POD , θ e ), X e POD (0) = X e POD + ζ X e POD , θ e = θ e + ζ θ e ,
where A e POD is computed from the finite element matrices M e , K e i and the vectors F e , G are projected on the POD space. Then, after discretization, this (small) system can be estimated by EKF or UKF (Appendix A). In doing so, the estimator is known to be stable, but this approach has a drawback: if the POD basis is not rich enough to capture the relevant phenomena, the reduced order dynamics may be a poor approximation of the full order one. Thus we may have a consistency problem in our way of modeling and approximating the system of interest. That is why we prefer another strategy which consists in solving the full order dynamics, and then applying the filter to the projection of the state vector on the POD basis. In other words, we apply RoUKF algorithm (Appendix B) with a reduced variable made of the parameters and the components of the augmented state vector on the POD basis:

X e r = α θ ∈ R r = R NPOD+d such that X e = X e ⊥ X e r (30)
Contrary to the first approach, the full model is preserved. The drawback is that it cannot be proved that the part of the error which is not filtered will not grow and then pollute the approximation. Thus, while the first approach could suffer from a consistency problem, the second one can suffer from a stability problem. Nevertheless, in all our test cases, this second strategy proved to be robust and accurate. By combining this strategy with UKF, the following algorithm is obtained to filter the electrical state. The initial condition projector is given by

L(0) = L ⊥ (0) L r (0) =   L ⊥ (0) L α (0) L θ (0)   =   0 1 1   .
Then we consider an adequate UKF sampling rule composed of weights and particles (see [START_REF] Julier | A new extension of the Kalman filter to nonlinear systems[END_REF] or Appendix A), we store the associated weights (α i ) in the diagonal matrix D α and precompute specific particles, here the so-called unitary sigma-points (i.e. with zero mean and unit covariance). More specifically, we consider the p = N POD + d + 1 unitary simplex sigma points, where N POD +d is the reduced space dimension [START_REF] Julier | Reduced Sigma Point Filters for the Propagation of Means and Covariances Through Nonlinear Transformations[END_REF][START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF]. We denote them by (I i ) 1≤i≤p and perform at each time step: a first sampling step which generates the particles identified by a subscript [i] ; a prediction step denoted by an additional superscript -and a correction step denoted instead by an additional superscript + .

1. Sampling:

             C n = (U n ) -1 X[i]+ ⊥n = X+ ⊥n + L X ⊥ n • C n • I i , 1 ≤ i ≤ p α[i]+ n = α+ n + L α n • C n • I i , 1 ≤ i ≤ p θ[i]+ n = θ+ n + L θ n • C n • I i , 1 ≤ i ≤ p (31a)
2. Prediction:

                         X[i]- n+1 = A n+1|n ( X[i]+ n ), 1 ≤ i ≤ p α[i]- n+1 = Π X[i]- n+1 , 1 ≤ i ≤ p X[i]- ⊥n+1 = X[i]- n+1 -Πα [i]- n+1 , 1 ≤ i ≤ p X- ⊥n+1 = p i=1 α i X[i]- ⊥n+1 α- n+1 = p i=1 α i α[i]- n+1 θ- n+1 = p i=1 α i θ[i]- n+1 (31b)
3. Correction:

                                             L X ⊥ n+1 = [ X[ * ]- ⊥n+1 ]D α [I [ * ] ] L α n+1 = [α [ * ]- n+1 ]D α [I [ * ] ] L θ n+1 = [ θ[ * ]- n+1 ]D α [I [ * ] ] Z [i]- n+1 = H n+1 ( X[i]- n+1 ) Z - n+1 = p i=1 α i Z [i]- n+1 Γ n+1 = [Z [ * ]- n+1 ]D α [I [ * ] ] U n+1 = 1 + Γ n+1 W -1 n+1 Γ n+1 X+ ⊥n+1 = X- ⊥n+1 + L X ⊥ n+1 M n+1 Γ n+1 M n+1 (Z n+1 -Z - n+1 ) α+ n+1 = α- n+1 + L α n+1 U n+1 Γ n+1 M n+1 (Z n+1 -Z - n+1 ) θ+ n+1 = θ- n+1 + L θ n+1 U n+1 Γ n+1 M n+1 (Z n+1 -Z - n+1 ) (31c)
where [I [ * ] ] is the matrix concatenating the (I i ) vectors side by side, and similarly for other vectors [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF]. From the last three corrections defined in [START_REF] Krishnamurthy | Patient-specific models of cardiac biomechanics[END_REF], the state correction reads

X+ n+1 = X- n+1 + L X ⊥ n+1 + ΠL α n+1 U n+1 Γ n+1 M n+1 (Z n+1 -Z - n+1 ).

State estimation in mechanics: Luenberger observers

To filter the mechanical state, a Luenberger approach is used, following [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF]. Consider the first two equations of (9) with

∂ t y = v, in Ω H 0 ρ∂ t v -∇ • (T ) = 0, in Ω H 0 ( 32
)
where the observations are given by

z m = y |ω obs m . ( 33
)
We then define the estimator by

     ∂ t ŷ = v + γExt ω obs m (z m -ŷ|ω obs m ), in Ω H 0 ρ∂ t v -∇ • ( T ) = 0, in Ω H 0 same boundary conditions as (9) ( 34
)
where γ is a scalar gain and Ext ω obs m is an extension operator typically given for any displacement field d by

d ext = Ext ω obs m (d) :              ∇ • (A : ε(d ext )) = 0 in Ω H 0 \ω obs m d ext = d in ω obs m A : ε(d ext )) • n = k s d ext , on Γ n (A : ε(d ext )) • n = 0, on ∂Ω H 0 \Γ n ( 35
)
with A the elasticity tensor coming from the linearization of T around 0 or a given trajectory at a given time t. On a linearized system, it can be proved that the the state error x = ỹ ṽ = y-ŷ v-v between the estimator and the target tends to zero [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF]. A consistent space discretization is given by

   ˙ Y = ˆ V + γ Ext(Z m -H m (X)), M m ¨ Y + K m ( Y , V ) = N m (36)
which converges to the real trajectory. Note that the proof of convergence after discretization is in general difficult and should be proved here using the internal viscosity to control the spurious high frequencies introduced by the discretization. Finally in a state-space form, there exists a Luenberger filter

G m L such that Ẋm = A m (X m ) + G m L (Z m -H m (X m )) ( 37
)
which converges for any initial error to an observed trajectory. This Luenberger filter, defined in the context of passive non-linear mechanics, was also shown to be robust to the introduction of the active part in the cardiac mechanics [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-mri using a biomechanical heart model[END_REF][START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF]. In this context, the Luenberger filter can be considered as a filter on the reduced space of the displacement and velocity field whereas the internal variables associated with the heart contraction are already stable. Nevertheless, note that the filter robustness has been demonstrated numerically but remains to be proved theoretically.

Estimation of one-way coupled systems

In this section we present an original strategy to aggregate the already defined filters of each submodel -namely the electrophysiological model and the mechanical model. This strategy is then justified by an estimation error analysis. To simplify the presentation and analysis of our method, we first consider the following one way coupling system with perfectly known parameters. Indeed, we have already seen how a joint state and parameter estimation can be added "in a second stage" once the state estimation is proved to be effective. Hence in practice the parameter identification will be considered but here, we can focus on the state as the main difficult aspect of the analysis. We thus consider Ẋe = A e (X e ), X e (0

) = X e + ζ e Ẋm = A m (X m , X e ), X m (0) = X m + ζ m (38)
and we propose to combine the Luenberger state filter on the mechanical part and the reduced order filter of Section 5.1 on the electrical part. We also recall that the observation operator is composed of two concatenated operators H e (X e ) and H m (X m ) as defined in [START_REF] Chapelle | An energy-preserving muscle tissue model: formulation and compatible discretizations[END_REF].

We can consider a RoEKF without loss by generality at the time-continuous level. Concerning the derivatives, we denote d m and d e respectively the partial differential with respect to X m and X e respectively, whereas we keep d when there is no ambiguity or when the differential is totali.e. with respect to all variables. From the general reduced order formulation of the EKF presented in Section 4.3 we introduce this time

L = L e L m such that                          Ẋe = A e ( Xe ) + L e U -1 L dH M(Z -H( X)), Xe (0) = X e Ẋm = A m ( Xm , Xe ) + G m L (Z m -H m ( Xm )) + L m U -1 L dH M(Z -H( X)), Xm (0) = X m Le = (dA e ( Xe ))L e , L e (0) = 1 Lm = (d m A m ( Xm , Xm ) -G m L dH m ( Xm ))L m +(d e A m ( Xm , Xe ))L e L m (0) = 0 U = L (dH) M(dH)L, U(0) = U e . ( 39
)
Note that in the dynamics [START_REF] Mirvis | Detection and localization of multiple epicardial electrical generators by a two-dipole ranging technique[END_REF] the Luenberger filter only applies on the mechanical part using only the mechanical data. However, the optimal filter strategy allows to benefit from both the electrical and mechanical data to correct the electrophysiology dynamics. This correction is then reverberated to the mechanical dynamics which also reads

Ẋm = A m ( Xm , Xe ) + G m L (Z m -H m ( Xm )) + L m (L e ) -1 ( Ẋe -A e ( Xe )),
when assuming L e invertible. Note that in the particular case where A e is linear the last expression simplifies into

Ẋm = A m ( Xm , Xe ) + G m L (Z m -H m ( Xm )) + L m • (L e ) -1 Xe ,

Convergence analysis

To ensure the convergence of the observer, L e is assumed to be invertible. Note that L e follows the tangent dynamics of the electrophysiological around the estimated trajectory of the model, starting from the initial condition L e (0) = 1. Therefore, even for a dissipative system, we can assume L e to be invertible but potentially ill-conditioned. Moreover the invertibility of L e can also be verified numerically in the specific case of concern. We recall that we denote by X = X † -X the estimation error. The variable X follows a non-linear dynamics. Then by linearization of this dynamics around the target trajectory we get a linear dynamics satisfied by linearized estimation error denoted by δ X = δ Xe δ Xm . In fact we have

       δ Xe = -U -1 L (dH) M(dH)δ X, δ Xe (0) = ζ e δ Xm = (d m A -G m L d m H m )δ Xm + (d e A)δ Xe + L m (L e ) -1 δ Xe , δ Xm (0) = ζ m . (40)
Here we introduce the change of variables

(δ Xm , δ Xe ) → (δη, δ Xe ) = (δ Xm -L m (L e ) -1 δ Xe , (L e ) -1 δ Xe ),
and obtain

       δη = (d m A -G m L d m H m )δη, δη(0) = ζ m • (L e ) -1 δ Xe = -U -1 L (dH) M(dH)L(L e ) -1 δ Xe -U -1 L (d m H m ) M m (d m H m )δη, δ Xe (0) = ζ e . (41)
The first equation corresponds to the dynamics of the linearized error studied for the mechanical system. Hence it converges to 0. Therefore the second term in the second equation tends to 0. The homogeneous part of the second equation can then be proved to converge to 0 if the following observability (42) condition is satisfied with our linear observation operator -namely, dH = H in our particular example. Namely we expect for all initial error δ Xe (0) that

∃(C, T ), T 0 H(L(L e ) -1 δ Xe ) 2 M ≥ C δ Xe (0) 2 U e , (42)
which can be at least verified numerically. In the last observability condition we see that L(L e ) -1 δ Xe represents the effect of a variation on δ Xe on both the electrophysiology and the mechanics. This effect is then observed through H. The observability is thus expected to be improved with respect to the situation where only the electrophysiology is considered. This will be confirmed numerically.

Practical algorithm

The correction steps of algorithms (B.5) (for mechanics) and ( 31) (for electrophysiology) are independent of the problem, and can thus be implemented by adopting again a master-slave strategy. In our electromechanical simulator, this yields a two-level interweaved master-slave strategy, the master handling the electromechanical problem becoming a slave of the sequential estimator algorithm. The tasks of the electromechanical master are the data exchange and the correction of unknowns and parameters. The estimation algorithm is summarized in Figure 4 for a generic time step. The first data transmission appearing in Figure 4 is done in order to update the initial value for each solver. This is handled by the electromechanical master.

Numerical illustrations

In this section, the methodology presented in Section 5.1 is applied to the electrophysiological problem of Section 2.1. The purpose is twofold: first, assess the reduced order Master RoUKF (Sampling) state filtering based on POD; second, illustrate the interest of estimating simultaneaously the state and the parameters.

[i] Particle Solver [p] Particle Solver [1] Particle Solver X[i]+ n , θ[i]+ n X[p]+ n , θ[p]+ n X[1]+ n , θ[1]+ n MASTER RoUKF (Correction) X[i]- n+1 , θ[i]- n+1 X[p]- n+1 , θ[p]- n+1 X[1]- n+1 , θ[1]- n+1
All the test cases of this Section are performed on the geometry of Figure 2, including the fibers in the conductivity tensors, and with the parameters typically used to generate healthy ECG [START_REF] Boulakia | Mathematical Modeling of Electrocardiograms: A Numerical Study[END_REF]. Synthetic data are generated by applying in (4) an external current I app during 25 ms. Two sets of measurements are considered: (1) 12-lead ECG only or (2) 12-lead ECG enriched with the extra-cellular potential at 8 epicardium points. The first set of measurements is of course the easiest to obtain. The rationale for the second set is to foresee what could be the benefit of including the measurements available from cardiac stimulators or implantable defibrillators. For both cases, these measurements are perturbed by an additive Gaussian white noise of standard deviation of 0.25mV (Figure 9). Note that another option could be to include endocardial extracellular potentials, if a catheter can provide them. We do not anticipate significant differences with respect to the results obtained with epicardial measurements.

ECG based state estimation of the electrophysiological model

Before addressing parameter estimations, a first test is run to assess the reduced order state filter, assuming that all parameters are perfectly known -for instance with τ in = 0.8 and τ out = 18. In other words, we study the problem of the state estimation alone applied only tp the electrophysiological model. The POD Basis used in the reduced order state filter [START_REF] Krishnamurthy | Patient-specific models of cardiac biomechanics[END_REF] was obtained, in this case, from snapshots generated with these parameters, among others. We choose to run the estimation simulation without any external current and compare it with the target simulation starting from t = 40 ms. When the estimation starts, the state variables are therefore significantly perturbed with respect to the direct simulation used to generate the synthetic data. In doing so, we introduce a significant state error in our estimation problem. Figure 5 compares the estimated mean value of the transmembrane potential v m and the ionic variable w with the target values, while Figure 6 compares the spatial distributions over the whole heart domain. We see that, despite the very large initial error, the estimator succeeds quite well in compensating for the lack of information about the initial stimulation and retrieving in time the accurate electrical state. Indeed, it is remarkable that with a correction only coming from the reduced state filter and applied from t = 40ms, the resulting estimated transmembrane potential and ionic current are in good agreement with their target values. Eventually, we point out that, as expected, adding a few measurements on the myocardium slightly improves the result. However this type of measurement clearly requires an invasive procedure [START_REF] Eldar | Transcutaneous multielectrode basket catheter for endocardial mapping and ablation of ventricular tachycardia in the pig[END_REF] problems is τ in = 1.5, τ out = 11. We point out that neither the values used to generate the synthetic data nor the initial ones are included in set of solutions used to generate the POD basis. In fact, the POD basis used in the reduced filter [START_REF] Krishnamurthy | Patient-specific models of cardiac biomechanics[END_REF] is computed from snapshots obtained with 4 sets of parameters τ in and τ out -presented in table 1 -which are in a neighbourhood of the expected value, but do not coincide with it as theoretically justified in [START_REF] Chapelle | A Galerkin strategy with Proper Orthogonal Decomposition for parameter-dependent problems -Analysis, assessments and applications to parameter estimation[END_REF]. Before launching the estimation, we proceed to a reparametrization of the form τ in = 2 θ1 τ † in and τ out = 2 θ2 τ † out , where τ † in and τ † out are given and θ 1 and θ 2 are the new values to be estimated. This is motivated by the fact that the parameters should be maintained positive during the estimation and the uncertainty variance is more naturally centered with respect to a power of 2 of initial parameter. We then compare two strategies: (1) parameter estimation only; (2) joint state-parameter estimation. In strategy [START_REF] Auroux | Back and forth nudging algorithm for data assimilation problems[END_REF], the component X r to be filtered is limited to the parameters θ, whereas in strategy (2), X r includes the parameters θ and the expansion coefficient α on the POD basis (see [START_REF] Kalman | New results in linear filtering and prediction theory[END_REF] The evolution of τ in and τ out is reported in Figure 7, and summarized in Table 2. In particular we clearly see in Figure 7 the convergence of the estimator on the parameters variable. As an additional illustration of the performance of the joint state-parameter estimator, we plot in Figure 8 8 leads of the corresponding standard body surface ECGs obtained with the initial guess of the parameters (curve named "ECG built from wrong initial guess"), with the parameters obtained after a simple parameter estimation (curve named "ECG after identification only"), with the estimated parameters obtained after a joint state-parameter estimation (curve named "ECG after joint state and param. estim."). All these curves should be compared to the one named "Target ECG" which corresponds to the measurements used for the estimation. Two comments are in order. First, we note that the ECG obtained with the "wrong initial guess" is indeed very different from the target. This means that our initial guess was far from the solution we looked for. Second, we see that the joint state-parameter estimation clearly outperforms the simple parameter estimation (as noted in a different context in [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF]). This means that only a joint state-parameter estimation allows us to produce an electrical state that is really compatible with the observations starting from realistic initial covariance -namely without forcing a priori the estimation to fit the data.

Full electromechanical data assimilation

In this section, we demonstrate numerically the efficiency of our complete state and parameter estimation chain for coupled models and, more importantly, we show that the estimation of a piecewise constant electrical parameter can be improved by enriching the electrical observations with kinematical observations. The parameter of interest is here τ close , which controls the plateau duration of the action potential of the cardiac cells and is supposed to be heterogeneous in the myocardium. In fact, τ close is assumed to take different constant values in 4 regions: inner part of the myocardium in the left ventricle, that will be called "endocardium" for simplicity; outer part of the myocardium in the left ventricle, that will be called "epicardium"; a thin region between the endocardium and the epicardium in the left ventricle, that will be called "M-cell"; and finally the right ventricle. Two kinds of synthetic observations are generated from a direct simulation: the electrocardiograms and the displacements of the myocardium. Hence for this numerical illustration, we have at our disposal a 3D field of displacement as it could ideally be processed from a 3D tagged-MRI sequence [START_REF] Rutz | Accelerated whole-heart 3D CSPAMM for myocardial motion quantification[END_REF]. However this is clearly an ideal situation from the mechanical point of view. Nevertheless, we believe that this illustration offers a clear insight into the maximum of information that could be earned with mechanical observations from the point of view of the electrophysiology model. These quantities are perturbed by an additive 10% gaussian white noise in time and space, proportional to the mean amplitude of the signal -see the corresponding signals in Figure 9. This type of noise is clearly illustrative as it is reasonable to choose independent time observations whereas the spatial distribution of the noise could certainly be more complex than a gaussian noise. In fact, real observations can even be spatially biased. However, this choice of noise is well adapted to evaluate a minima the robustness to noise of an estimation method. Facing these noises we follow the recommendation made in [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-mri using a biomechanical heart model[END_REF] Figure 10 shows the time evolution of τ close during the sequential estimation in each region for electrical measurements only and for the electromechanical -ECG plus displacements -available measurements. The benefit of the electromechanical measurements is striking, especially in the M-cell region. Then we present the final estimated values and the error with respect to the target values in Table 4. We also report the final estimated values for electromechanical measurements after three heart beats where the estimation error is even more reduced. With the electrical measurements only, we were not able to run three heart beats, because the estimated parameters diverged too much from physiological values. On the contrary, with electromechanical measurements, the results keep improving along the three heart beats. Thus, it clearly appears that the estimation is much more accurate and robust when the electromechanical measurements are taken into account. Our data assimilation procedure results can also be evaluated in the light of ECGs that the model can produce. Indeed, we plot in Figure 12 the reference ECG (in dashed black), the ECG corresponding to the initial guess (in black), the ECG obtained from the estimated values with the electrical measurement (in gray) and finally the ECG obtained from the estimated values with the electromechanical measurements (in dashed gray). We clearly see that the ECG corresponding to the parameter initial guess has a completely false T-wave, whereas after the estimation procedure the ECGs are much closer to the reference, especially with electromechanical measurements.

In order to quantitatively understand the relevance of using the mechanical observa- tions to estimate an electrical parameter, a sensitivity analysis is performed. To simplify the presentation of this last calculus, let us consider a time continuous ROEKF for the parameter estimation. We recall that the sensitivity of the model with respect to its parameters is directly estimated with L X in the ROEKF filter [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF]. In our practical case we recall that L X is a matrix which consists of 4 columns vectors (one for each region of the Table 3) of dimension the number of the electrical model degrees of freedom plus the mechanical model degrees of freedom. Therefore, the measurements sensitivity with respect to the parameters is estimated by HL X . Note that we have already seen such sensitivity terms in the observability condition [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF] introduced in the convergence study of our coupled estimator. In our practical case, we focus on the sensitivity s e , resp. s m , of the first lead of the ECG, resp. the mean value of the measured displacements norm in the myocardium, with respect to τ close , and we consider normalized quantities. We denote by (Z e ) 1 the first index of Z e where the first lead of the ECG is gathered. We In practice when using the ROUKF the computations are very similar but time integrations is replaced by time iteration summations whereas HL(t) are replaced by Γ n+1 (L θ n+1 ) -1 computed from [START_REF] Krishnamurthy | Patient-specific models of cardiac biomechanics[END_REF]. Our two normalized sensitivities are represented in Figure 11 for the four different regions. In the epicardium and in the right ventricle, the sensitivity of the electrical measurement is higher than the mechanical one, as expected. But interestingly, the electrical sensitivity is lower than the mechanical one in the M-cell and endocardium regions. This is particularly the case here where we observe the displacement in the whole heart. However recent advances in tagged-MRI allows to expect that intra-myocardial kinematics will be available in the future whereas it is out of reach for non-invasive electrical measurement. In addition, the mechanical observation is affected by τ close over a longer time window. Hence, the mechanical measurements therefore increase the time interval during which the observer has the opportunity to correct the value of τ close .

In conclusion, this last example is a rich illustration of our sequential strategy, with a state filter based on a Luenberger approach for the mechanics and on a reduced order optimal filter for the electrophysiology and the parameters. In particular it allows us to

V6

ECG built from wrong initial guess ECG after identification with ECG obs. only ECG after identification with electromechanical obs. Target ECG Figure 12: ECGs obtained before and after estimation, compared to the ECG used as a measurement numerically justify the importance of having a correction for every type of uncertainties, especially for state uncertainties. Moreover we numerically demonstrate the interest of taking advantage of the multi-physics nature of a problem to identify some parameters -here τ close accounting for the plateau duration -which have an impact on the coupled physics.

Conclusion

In this work, a complete framework for the joint state and parameter estimation of a complex multi-physics problem has been presented. The strategy consists of Luenberger observers and reduced order optimal observers. The robustness of the approach is ensured by the fact that each component of this system is corrected by the data. In particular, the importance of state estimation when performing a parameter identification has been illustrated. For weakly coupled systems we have justified that our filter combination and aggregation lead to a reduction of the estimation error. It is expected that the same type of results could be proved for fully coupled problems in the future. This general strategy has been applied to an electromechanical model of the heart but other coupled systems of interest in various fields of application can be estimated with this approach. In this framework, we have also shown that the use of mechanical data could considerably improve the observability of the problem.

The direct electromechanical model is clearly more demanding than the direct electrophysiology model. However, we have shown that the associated inverse problem which consists in identifying its parameters from all the available data -mechanical and electrical -is better defined than the pure electrophysiology inverse problem.

In conclusion, this work provides a promising new strategy to address the classical inverse problem of electrocardiography, and suggests a new way to reduce its well-known ill-posedness. The numerical results, based on synthetic data, are presented primarily to validate the methodology. This validation being successfully achieved, an important issue to be addressed in future works concerns the adequacy and accuracy of the models. Another important perspective is the validation against real data. Using real data will necessitate a complete patient workflow anatomical ECG and Tagged-MRI (or at least Cine-MRI based on the results of [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-mri using a biomechanical heart model[END_REF]) are acquired and post-processed. In this respect, we point out that our strategy is indeed quite demanding in term of data acquisition and processing. This must be ultimately compared -in term of estimation robustness and accuracy versus difficulty-to existing methods such as electrocardiographic imaging (ECGI) [START_REF] Ramanathan | Noninvasive electrocardiographic imaging for cardiac electrophysiology and arrhythmia[END_REF] where only ECG-like measurements are necessary with special vests recording up to 224 body-surface potentials.

where [I [ * ] ] is the matrix concatenating the (I i ) vectors side by side, and similarly for other vectors.

The analysis of this estimator can also be found in [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF] and relies again on the fact that the linearized error satisfies exactly the same dynamics as the linearized error of the RoEKF. Therefore the demonstration by linearization is directly obtained. The combined practical simplicity and efficiency of such parameter estimator in comparison to other adaptive observers have made this approach popular for real case parameter estimation problems [START_REF] Bertoglio | Sequential parameter estimation for fluid-structure problems. application to hemodynamics[END_REF][START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-mri using a biomechanical heart model[END_REF][START_REF] Marchesseau | Personalization of a Cardiac Electromechanical Model using Reduced Order Unscented Kalman Filtering from Regional Volumes[END_REF][START_REF] Xi | Myocardial transversely isotropic material parameter estimation from in-silico measurements based on a reduced-order unscented Kalman filter[END_REF].

Figure 1 :

 1 Figure 1: Master -Slave coupling

Figure 2 :

 2 Figure 2: Mechanical mesh (Left, 4907 nodes, 18193 elements) and electrical mesh (center, 108112 nodes, 541994 elements) and thorax mesh (right, 229782 nodes, 1250072 elements)

Figure 3 :

 3 Figure 3: Direct simulation results in long axis view -deformed mesh with interpolated transmembrane potential from electrical refined mesh.

Figure 4 :

 4 Figure 4: Estimation algorithm: "particle solver" denotes the electromechanical solver presented in Section 2.3

Figure 5 :

 5 Figure5: Estimation of the transmembrane potential vm (Left) and the ionic variable w (Right). The direct simulation was performed with an external current applied during 25ms. The curves represent space-averaged quantities.

Figure 6 :

 6 Figure 6: Space distribution of the transmembrane potential for the direct simulation (top) and the state estimation with ECG measurements (bottom). There is no initial activation in the bottom simulation in order to generate a strong error in the initial condition. The correction of the filter is applied from t = 40 ms. Note the good agreement of both simulation from t = 55ms.

Figure 7 :

 7 Figure 7: Comparison on τ in (Left) and τout (Right) estimation, with (black line) and without (gray line) reduced state filtering. These curves clearly show that reducing the uncertainties on the state actually improves the identification of parameters

Figure 8 :

 8 Figure 8: ECGs after state and/or parameter estimation (zoom on a time window of 400 ms)

Figure 9 :

 9 Figure 9: I lead electrocardiogram with additive Gaussian noise (left) and displacement of a mesh point with additive Gaussian noise (right)

τ

  with real data to choose a time-discretized observation norm built from the inverse of the standard deviation δ obs by ∀n, M e n = ∆t e (δ e obs ) -2 1, and , M m n = ∆t m (δ m obs ) -2 1. with δ e obs 10% * 2.5 mV = 0.25mV and δ m obs 10% * 10 mm = 1 mm. As in Section 6.2, the unknowns are reparametrized as τ close = 2 ϑ τ † close , where τ † close corresponds to the initial guess. The values to be estimated and the initial guess are gathered inTable 3. close RV estimation Mean value estimation from ECG Standard deviation estimation from ECG Target Value Mean value estimation ECG+displ. Standard deviation estimation from ECG+displ.

Figure 10 :Table 3 :Figure 11 :

 10311 Figure 10: τ close estimation on endocardium, m-cell, epicardium and right-ventricle regions, comparison ECG observations vs ECG+mechanical observations

thus define s e (t) = τ close ze 1 (

 1 H e L e (t)) 1 , with ze 1 e (t)) 1 | dt, and T is typically the heart beat duration. Identically we defines m (t) = τ close zm ((H m L m (t)) M m H m L m (t))

  .

	Slave Elec. Solver		Master Elec.-Mech.		Slave Mech. Solver
	Elec. Solver t n	Start Iter.	Initialization	Start Iter.	Mech. Solver t n
	•	∆t e
	•			n = n + 1	∆t m
	•
	t n + ( ∆tm ∆te )∆t e t n+1 = Elec. Solver	v n+1 m	Interpolation	v n+1 m	t n+1 = t n +∆t m Mech. Solver

Table 1 :

 1 . Set of parameters used to build the POD basis.

	Simulation	1	2	3	4
	τ in 0.8	0.8	1.2	1.2
	τ out 18.0 12.0 18.0 12.0

  ).

	Param. Target	A priori	Only Param.	Joint State &
		Value	(%Error)	Identification Parameter Estim.
	τ in	1.0	1.5(50%)	2.19(118.8%)	1.02(1.99%)
	τ out	16.0	11.0(31.25%) 6.96(56.48%)	15.15(5.3%)

Table 2 :

 2 Target vs A priori values for τ in and τout

Table 4 :

 4 Identification of τ close in the various regions and for the different scenarios

Preprint submitted to ElsevierDecember 6, 2014

the observer is called optimal due to the fact that it is associated with an optimal control problem.

Acknowledgement

The authors would like to thank the anonymous reviewers for their valuable comments and suggestions.

Appendix A. Time-discrete EKF and UKF filters

The time-discretization of the optimal estimator presented in Section 4.2 is based on the principle that the optimality should be conserved at the discrete level. Let us first denote a stable and consistent discretization of the original model by

and define a discrete-time functional

with M k = ∆tM when a fixed time-step is considered. This choice of discrete observation norm M k ensures that the discrete-time functional J N is consistent with respect to the continuous-time functional J . The discrete-time optimal filter is then defined by Xn = X n,[argmin Jn] . The discrete-time counterpart of [START_REF] Gulrajani | Moving dipole inverse ecg and eeg solutions[END_REF] for the optimal linear estimator can be deduced with a prediction-correction time scheme. We denote by a superscriptthe prediction computations and by a superscript + the correction computations. 1. Prediction:

This algorithm can be interpreted as prediction-correction time discretization of [START_REF] Gulrajani | Moving dipole inverse ecg and eeg solutions[END_REF].

Another way to propose an efficient estimator is to replace at the time-discrete level the tangent operators by a finite difference interpolation based on the computation of the original operator on numerous sampling points. This is the case for Ensemble Kalman Filter (EnKF) [START_REF] Evensen | Data Assimilation: The Ensemble Kalman Filter[END_REF] or the Unscented Kalman Filter (UKF) [START_REF] Julier | A new extension of the Kalman filter to nonlinear systems[END_REF]. In this article we focus on the second one which is a discrete-time estimator based on sampling particles called sigma-points helping to replace the tangent computations. Let us introduce the so-called unitary sampling points I [i] and weight α i with the following rules

so that, at each time step, the sigma-points can be generated around the estimated values based on the covariance estimation by

then we compute the prediction

and the corresponding observations

The gain is then defined by

so that, with the covariance

This filter has the great advantage to be much easier to compute numerically since, contrary to EKF, the tangent operators computations are no longer required.

Appendix B. Reduced-Order filters: time discrete RoEKF and RoUKF

The reduced-order filtering concept presented in Section 4.3 can be directly applied to the time and space discretized versions of the equations, leading to a time-discrete RoEKF:

1. Prediction:

2. Correction:

This algorithm has been applied in [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF] for parameter identification by reducing the uncertainty space to the parameter space. Indeed, the extension L is initially decomposed into

Using the fact that the parameters dynamics is null, we can easily proved that ∀t > 0, L θ = 1.

Then, by direct computation, we have the continuous-time formulation of the RoEKF in a parameter identification context

In (B.3), we see that the filter correction appears in the parameters dynamics. Then, the correction on the parameters is reverberated to the state by the mean of L X which can be interpreted from its proper dynamics as the sensitivity of the model with respect to the parameters. Therefore even for a reduced order strategy on the parameters the global filter corrects the two components of the model, namely the parameters but also the state.

The same principles can be applied after the time discretization of the model to get 1. Prediction:

starting from the same initial conditions. We have presented the Order EKF and have recalled its convergence. Now a legitimate question is to know if a reduced order strategy can be applied to other types of approximation of the optimal filtering approach. The answer was given for the UKF filter in [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF] where a reduced order version (RoUKF) was derived. Moreover it was also applied to a specific case where the initial uncertainty is reduced to the parametric space. In the general contexti.e. without particularizing the state and parameter dependency -the RoUKF can be formulated as follow.

Given an adequate sampling rule, we store the corresponding weights (α i ) in the diagonal matrix D α and precompute specific unitary simplex sigma-points (I i ) 1≤i≤p (i.e. with zero mean and unit covariance) with p = d + 1 since the reduced space corresponds here to the parametric space. The algorithm consists of the following three steps computed recursively:

1. Sampling:

3. Correction: