Sequential data assimilation for mechanical systems with complex image data: application to tagged-MRI in cardiac mechanics

Tagged Magnetic Resonance images (tagged-MRI) are generally considered to be the gold standard of medical imaging in cardiology. By imaging spatially-modulated magnetizations of the deforming tissue, indeed, this modality enables an assessment of intra-myocardial deformations over the heart cycle. The objective of the present work is to incorporate the most valuable information contained in tagged-MRI in a data assimilation framework, in order to perform joint state-parameter estimation for a complete biomechanical model of the heart. This type of estimation is the second major step, after initial anatomical personalization, for obtaining a genuinely patient-specific model that integrates the individual characteristics of the patient, an essential prerequisite for benefitting from the model predictive capabilities. Here, we focus our attention on proposing adequate means of quantitatively comparing the cardiac model with various types of data that can be extracted from tagged-MRI after an initial image processing step, namely, 3D displacements fields, deforming tag planes or grids, or apparent 2D displacements. This quantitative comparison—called discrepancy measure—is then used to feed a sequential data assimilation procedure. In the state estimation stage of this procedure, we also propose a new algorithm based on the prediction–correction paradigm, which provides increased flexibility and effectiveness in the solution process. The complete estimation chain is eventually assessed with synthetic data, produced by running a realistic model simulation representing an infarcted heart characterized by increased stiffness and reduced contractility in a given region of the myocardium. From this simulation we extract the 3D displacements, tag planes and grids, and apparent 2D displacements, and we assess the estimation with each corresponding discrepancy measure. We demonstrate that—via regional estimation of the above parameters—the data assimilation procedure allows to quantitatively estimate the biophysical parameters with good accuracy, thus simultaneously providing the location of the infarct and characterizing its seriousness. This shows great potential for combining a biomechanical heart model with tagged-MRI in order to extract valuable new indices in clinical diagnosis.


Introduction
Cardiac biomechanical modeling has made tremendous progress over the past decades, and some accurate models are now available to represent the complex deformations of the organ-among other quantities of interest-over full heartbeats, frequently based on multi-physics and multi-scale formulations, see e.g. [1,2] and references therein.
As for all natural systems-as e.g. in geophysics-a great challenge consists in dealing with the many unknown or uncertain quantities-initial conditions, boundary conditions, and various physical parameters-that need to be prescribed for running model simulations [3]. In this work, we decide to rely on a data assimilation strategy [4] to estimate the uncertain quantities while allowing predictive simulations.
Concerning the specific problem of estimation in cardiac biomechanical modeling, difficulties arise from both (1) the complexity of the models considered, and (2) the nature of the available measurements, often relying on medical imaging [5]. An effective estimation methodology has been proposed by [6] for this type of model, based on a so-called sequential approach-also known as observer method. In this approach, the dynamical model is corrected at each time using the computed discrepancy between the current simulation and the actual measurements, see also [7]. This strategy was designed to be applicable to measurements concerning displacements, whether they be given internallyin a sub-region of the system-or on a boundary or a part thereof. It was also shown to be extendable to data consisting of segmented surfaces as obtained by processing various types of medical imaging dynamical sequences.
In this paper, we focus on estimation based on data provided by tagged-MR imaging sequences [8,9]. Tagged-MR is generally considered to be the "gold standard" in cardiac imaging, in particular as regards the assessment of so-called "cardiac mechanical indicators", namely, indicators pertaining to displacements, strains, and volumes [10]. As a matter of fact, tagged-MR images visualize the deformations of grids associated with the actual tissues, which is of course most valuable for clinical purposes, both from a qualitative standpoint as assessed by the physician's eye, and with a view to obtaining such quantitative indicators. However, the problem of extracting actual 3D material displacements from a tagged-MR sequence gives rise to serious difficulties [11,12]. In fact, in many cases only 2D "apparent" displacements are obtained, which introduces specific errors in the displacement-based quantitative indicators, in addition to usual inaccuracies pertaining to image processing. Of course, these difficulties are also of concern when extracted displacements are to be used in an estimation setup, hence this justifies looking more closely into tagged-MR modalities to devise and analyze strategies to adequately employ them for estimation purposes. In this regard, the contributions presented in this communication are twofold.
First, we propose a systematic approach to incorporate within an estimation framework a wide range of data, potentially obtained from prior processing steps applied on tagged-MR images. These data vary in their nature, covering the cases of: full mechanical displacements in a subdomain; sequences of deforming tag planes and tag grids; and 2D apparent displacements. Extracting state-and parameter-corrections from these data is an intricate task. To address this challenge, we devise for every case relevant discrepancy measures. The soundness of our approach is corroborated by a complete mathematical analysis of the state observer in an idealized fully linear case, provided as complementary material in Appendix A.
Secondly, we propose a relevant time-discretization scheme for the state observer, which is a particularly crucial point in the context of sequential data estimation. This scheme is built upon a "prediction-correction" strategy, where the former step corresponds to genuine model time marching, and the latter to discrepancy-based adjustments. This clean decomposition of these steps offers numerous practical benefits. Additionally, we are able to provide evidence that the obtained time-discrete observer retains the convergence properties of the time-continuous observer, with rigorous proofs detailed in Appendix B.
The outline of the paper is as follows. In the forthcoming section we recall the main principles underlying the design of observers, and we provide a quick overview of the mechanical model of a beating heart, as an example of a model formulation. The next section is dedicated to describing the potential information extracted from tagged-MR images and to proposing-for each type of data-the discrepancy measures. We then address the issue of space and time discretization of the observer in order to perform joint state and parameter estimation. Finally, in the last section we present several numerical experiments in which we performed parameter estimation based on synthetic measurements.

Principles of sequential estimator design
The aim of a sequential estimator-also called observer-is to approximate a real trajectory, in spite of various uncertainties, using the knowledge provided by the measurements obtained on this specific real trajectory. Let us consider a real trajectory y ref (t), t ∈ [0, +∞), belonging to the so-called state space Y and solution, in our case, of a-possibly infinite-dimensional-dynamical system summarized in the state space forṁ with an uncertain initial condition y ref (0) = y 0 + ζ y , where y 0 is a known a priori and ζ y is the uncertain part in the initial condition. Therefore, any simulation of y-based on the discretization of the dynamical system-starting only from y 0 will be affected by the propagation of this error made in the initial condition. To circumvent this difficulty, we can benefit from the measurements at our disposal on the trajectory. We denote by z these measurements-also called observations and belonging to the observation space Z-which are assumed to be generated by a mapping H on the real trajectory, up to additional measurement errors z = H(y ref , t) + χ.
The observer denoted by y is a system that starts from the only part known in the initial condition-namely y 0 -and uses in time the available measurements z to generate a trajectory y(t), t ∈ [0, +∞) that converges to y ref as fast as possible. Therefore, simulating y instead of y from y 0 gives a better approximation of the targeted system.
The main categories of observers addressed here are computed by a feedback law based on the measurements in the form ⎧ ⎨ ⎩˙ y = F( y, t) + G(z − H( y, t)) where G is called the gain operator, also referred to as filter. The dynamics of y is corrected when a discrepancy is observed between the actual measurements z and the measurement H( y) that would have been produced by y. This discrepancy J( y, t) = z − H( y) is also called innovation since it not only expresses an observation error, but also a source of improvement for the observer. We point out that with certain types of measurementsas is typically the case with image-based observations-it is sometimes difficult to define an adequate observation operator but easier to directly compute a discrepancy [6]. This is not a problem for the observer definition since only the discrepancy appears. In a fully linear situation, namely, when the dynamics is linear with F(y, t) = A(t)y + R and H(y, t) = H(t)y, the most well-known gain operator is given by the Kalman gain, see e.g. [13,14] and references therein. This operator is expressed as G(t) = P(t)H(t) * where P is an operator following the Riccati evolution equatioṅ P = AP + PA * − PH * HP, P(0) = P 0 , and H * is the adjoint of H. Although P is computable for any dynamics operator A, it leads after spatial discretization to a discrete operator which is intractable in practice. This phenomenon has been known for decades and called "curse of dimensionality" [14,15]. Therefore, for specific dynamics other types of gains have been investigated as initiated by [16]. They are based on the fact that, when computing the estimation error y = y ref − y, we get in a fully linear setting the following dynamics Hence, G should be designed to stabilize the estimation error dynamics operator A − GH, so that the homogeneous system tends to 0, namely, In the presence of noise in the measurements, this would also control the error dynamics. This strategy is referred to as the Luenberger observer or nudging-see [17] for a survey. For the elastodynamics system-a particular case of second-order hyperbolic systems- [6] has shown that a very simple choice of with γ a scalar coefficient can be sufficient. In a nonlinear configuration, fewer theoretical results are available. However, an accepted strategy is to replace in the gain the use of the adjoint of the observation operator, namely H * , by the adjoint of the tangent operator DH( y) * around the estimated trajectory. Therefore for small errors we can expect that the linearized error around the trajectory is stable.
One last fundamental aspect that we need to describe in this introduction to observer design is how parameter estimation-also called identification-can be carried out. Let us denote by θ the uncertain parameter to be identified. Note that θ may be a vector of components or even a distributed field. The main idea is to introduce an augmented state vector and dynamics operator such that we still have ẏ = F( y, t). Then, a Kalman observer can be directly defined on this augmented model leading to a covariance operator and gain However, it is more intricate to define a relevant Luenberger observer for the augmented system as the observations are frequently linked to the parameters through the state only. Therefore, there is little hope that γ H * will lead to an efficient gain. An alternative strategy was proposed by [18] as a generalization of the adaptive filtering strategy of [19,20]. The idea is to retain the Luenberger observer on the state while using a Kalman-like gain on the parameters. This strategy can be very effective in practice, since it is common to consider a parameter described much more coarsely than the state discretization, thus alleviating the curse of dimensionality associated with optimal filtering. The complete observer reads where U −1 is a reduced covariance operator on the parameter space and L is a "sensitivity" operator from the parameter space to the state space. We see in the dynamics (2) that the state gain is the combination of the Luenberger gain and a gain directly inferred from the parameter filter so that In [18] the convergence of the complete observer is also established-at least in a linear configuration-based on the idea that the Luenberger state observer reduces the uncertainty to the parameter space where the optimal filter operates. Moreover, the effectiveness of this approach has already been applied to biomechanical identification problems by [7,21].

Example of model formulation
We consider here a model of beating heart involving a large displacement solid formulation with active stresses driven by an electrical input. Let us now introduce some notations in order to summarize the model equations. We denote the heart domain by (t) at any time t. This domain is the image of a reference configuration 0 through the solid deformation mapping where u denotes the solid displacement, so that the solid velocity is given by v =u. The deformation gradient tensor F is given by Furthermore, we introduce the right Cauchy-Green deformation tensor C = F · F . We finally recall that the Green-Lagrange strain tensor denoted by e is defined by Regarding the mechanical quantities notation, we denote by ρ the tissue mass per unit volume, and by σ the Cauchy stress tensor associated with the deformed configuration. In the reference configuration, we respectively define the associated first and second Piola-Kirchhoff stress tensors as The constitutive law can be considered as a nonlinear rheological combination of a passive part and an active part T = p + a . The passive part p is described by a hyperelastic law of potential W and a viscous component chosen proportional to the strain rateė Concerning the hyperelastic law, there exists some experimental evidence-based on detailed ex-vivo tri-axial shear testing-in favor of a complete orthotropic passive behavior [22], with a so-called sheet structure providing a second privileged direction, namely, after the muscle fiber direction. However, the sheet direction cannot be characterized invivo for patient-specific modeling purposes. Moreover, various studies have shown good agreements of transversely isotropic models with experimental data obtained at the organ level, see e.g. [23,24]. We thus consider a transversely isotropic law of exponential type earlier proposed by [25], and inspired from the fully orthotropic model of [26], viz.
where J 1 is the standard first reduced invariant, J 4 is the reduced invariant accounting for the anisotropy of the material in the fiber direction τ, namely J 1 = J − 2 3 tr(C), J 4 = J − 2 3 τ · C · τ, and κ is the bulk coefficient. For the active part a , we rely on the model proposed in [27], with internal variables defining the active strain e c , the active stiffness k c and the associated active stress τ c , along the fiber direction τ in a chemically-controlled constitutive law describing myofibre mechanics [27,28]. Therefore, we have a = σ 1D (e c , k c , τ c )τ ⊗ τ. We finally end up with the following second Piola-Kirchhoff stress tensor Concerning the boundary conditions, following [7] we model the interactions with the surrounding organs by visco-elastic boundary conditions on a subpart of the boundary, which gives in the reference configuration T · n = k s u + c s v on n , where n denotes the surface normal in the reference configuration. Regarding the pressure load, we consider a uniform following pressure on the left and right endocardium easily written in the deformed configuration σ · n t = −p v,i n t on n,i (t), i = {1, 2}, where, here, n t denotes the normal of the deformed configuration boundary. Finally, the complete mechanical model Together with the internal variable dynamics given in [27], it constitutes a general definition of the dynamical operator denoted by F in our above summarized description for a state y corresponding to (u, v, e c , k c , τ c ).

Filtering data available from imaging modalities
For assessing the physiological condition of a patient's heart, physicians usually seek standard indicators such as mass, volume or ejection fraction. Additionally intra-myocardial deformations are of great importance to assess the cardiac function in a localized manner. Even though the former three global indicators can be obtained using various types of medical image modalities-e.g. cine-MRI-the latter are more intricate to capture by standard procedures. Magnetic resonance imaging with tissue marking-referred to as tagged-MRI-was introduced in the late 80's [8]. By non-invasively imprinting a pattern in the acquired images-through specific magnetization of the tissue-it reveals myocardial deformations. Various types of tagged image modalities have arisen since its inception. They differ by the orientation, the temporal persistence or even the shape of the pattern. For instance, the SPAcial Modulation of Magnetization (SPAMM) modality-introduced by [9]-generates a grid-like pattern, whereas the first tagging images included a radial pattern-see e.g. [29]. The temporal persistence of the pattern in SPAMM images covers the complete heart systole. Figure 1a shows an example of a SPAMM image (in short axis view) at marking time. We observe the regular pattern within the image domain. Figure 1b shows the same image slice obtained during contraction. In this second image we observe the deformation of the originally regular pattern subject to the material displacements.
Even though SPAMM images are the most popular type of tagged-MRI, other modalities exist. For example, we can cite the DANTE sequence-initially introduced by [30]-which a 2D slice at marking time b 2D slice during systole provides a thinner tag grid pattern. Another example is the CSPAMM modality that aims at decreasing the tag pattern fading using two sequences of SPAMM images-see [31].
Combinations of 2D images, e.g. by superposition in a orthogonal direction, are historically the first type of tagged image modalities that was proposed. However, they are naturally affected by through-the-plane motion. This particular aspect will also be referred to as "2D apparent displacement" in the following. It corresponds to the shift in the position of the intersection of the deforming material with the image plane. Note that this displacement is neither a full displacement measurement nor a projection of the displacement onto the image plane. To circumvent this limitation, later works have led to the production of complete three-dimensional tagged-MRI-see [11]. 3D tagging (3D SPAMM) is an imaging modality of major interest since it can provide truly three-dimensional information on the heart strain.
From the most direct type of data to more complex observations, our work is dedicated to the design of observers based on 2D or 3D SPAMM images. Such an endeavor requires-see "Principles of sequential estimator design" section-to be able to compute the discrepancies between the various types of pre-processed data and the model. First, we assume that this image processing step leads to the reconstruction of the fully three dimensional deformation of the heart-from 3D SPAMM for instance [11,32], or from the collection of various 2D SPAMM [33,34]. In this context, following [6] we propose an efficient way to assimilate this direct displacement measurement. However, this may introduce a new source of error pertaining to displacement tracking, in addition to those inherent to the imaging modality per se. Hence, we further consider three distinct situations aiming at gradually decreasing our demands on this prior processing step. To start with, we propose a discrepancy measure based on the assumption that we are able to reconstruct the tag planes fitting the tag pattern [35][36][37]. In the case of two-dimensional images, obtaining these surfaces may require a complex interpolation scheme in the image transverse direction. Therefore, in a second step, we consider the case of 2D tag grids lying within the image planes [38]. In a final step, we propose means of comparison between the model and 2D apparent displacements extracted from 2D MR images, see [39][40][41][42].

Filtering 3D displacements from 3D grids
By constructing the tagging pattern in three directions, 3D SPAMM is a powerful image modality that potentially leads to a reconstruction of the complete three dimensional heart motion [11]. For instance, [12] proposed to adapt the HARmonic Phase (HARP) method-which performs tag patterns tracking by analyzing the frequency contents of the image-to 3D images in order to extract these data. However, this modality suffers from long acquisition times, requiring multiple breath-holds from the patient. Note that recent works [43] have shown that this difficulty can be partially overcome, typically using signal under-sampling. Another technique to extract three dimensional displacements of the heart from SPAMM images is to acquire two orthogonal sets of 2D images in short and long axis-see for instance [33,44] or [45]. It should be noted that this process is likely to suffer from slice misregistration and through-the-plane motion.
However, as a first step, it seems reasonable to assume that the observations take the form of the 3D tissue displacements with a resolution corresponding to the tag pattern spacing. In a (forthcoming) discrete setting, this entails the use of a space interpolation operator between the mesh model and the observation region. Nevertheless, in a continuous formalism, we can assume that we have at our disposal some measurements of the displacements in a subdomain ω 0 of the geometry. We introduce the observation operator where Y u = H 1 ( 0 ) 3 is the displacement space. Note that H apply to the state space Y gathering the displacement space and the velocity space (and eventually additional variables in more general mechanical configuration). By contrast, H |Y u is the same operator of input space restricted to the displacement space Y u . We then need to specify the observation space Z. One possible choice is to consider L 2 (ω 0 ) 3 , the space of square-integrable vector fields. However, this space does not characterize the maximum amount of information we have on the system, since the observations typically comes from a displacement in H 1 ( 0 ) 3 , the space of square-integrable vector fields with square-integrable first derivatives. We should rather consider Z = H 1 (ω 0 ) 3 , and we propose in Appendix A a more complete mathematical justification of this choice. Hence, following [6], we introduce a convenient way to define a norm in this space. Let us consider the following extension where σ lin denotes the stress tensor given by linearized isotropic elasticity. In particular, we denote the corresponding linear constitutive law by where ε denotes the usual linearized strain tensor. The boundary conditions on ∂ 0 are adequately obtained from (3). In the following, we will denote the extension operator by Note that an equivalent variational characterization of the extension is given by where the energy dot-product is here defined by We can prove-see Appendix A for a similar dot product (·, ·) E 0 -that It is now possible to define the adjoint of the observation operator that is needed in (1). We find in [6] and Appendix A that H * is given by Note that in the rest of the document, by a slight abuse of notation, we will denote by H the operator applying either on the state space or on the displacement field extracted from the complete state y. We can now define, in a continuous formalism, the state observer. As recalled in "Position of the problem" section, this state observer is the first ingredient of our state-parameter estimation procedure as the state is corrected by physically-based feedback and the parameters by a Kalman feedback. Indeed, following [6], the state observer corresponding to (1) with G = γ H * is given in strong formulation by: In order to justify why such a simple feedback is in fact very effective for controlling state errors in our formulation we propose: (1) a complete mathematical analysis in a simplified elastodynamics configuration in Appendix A; (2) an energy estimate of the estimation error proving at least the decrease of the estimation error in a general framework. Indeed, let us compute in the simplified case of linearized visco-elasticity without activation internal variables, namely The estimation error satisfies the following weak formulation, for any v ∈ Y u , with the additional observer-based relatioṅ assuming zero measurement error to fix the ideas. Weighing the latter relation by u and using the energy dot-product yields where we have used the orthogonality property (5). We can now substitute this expression in the above variational formulation applied with the test function v = v, which gives d dt where v 2 K = (ρ v, v) L 2 ( 0 ) 3 denotes the kinetic energy of the error. We can see that the total energy of the error-namely, elastic energy of the deformation plus kinetic energy of the velocity-decreases. Due to the observer correction term, we can expect a faster stabilization than with the sole natural dissipation.

Filtering tagged-MR planes and grid
The limitations of 3D SPAMM that we have previously mentioned make this imaging modality difficult to use in clinical routine. In most clinical cases only 2D tagged-MRI is available. These datasets can be treated plane by plane to extract apparent displacements. Prior to proposing a corresponding observer for this type of data, we assume that a first step of image processing leads to the construction of geometrical objects taking the form of tag planes or tag grids and following in time the deformations of the tag patterns-see [36][37][38] for examples of tag planes constructions and [38,46,47] for tag grids.

Extension of surface data
For the definition of the spaces and norms associated with the discrepancy measures needed in this section, following [6] and the approach detailed in the previous section, we will use an extension operator mapping data provided on a surface to the whole solid domain. More precisely, we denote by S 0 a surface embedded in the reference domain 0 and e a vector field defined on S 0 , with (e 1 ⊥ , e 2 ⊥ ) defined so that (e 1 ⊥ , e 2 ⊥ , e) gives an orthonormal basis at any point in S 0 . We define the extension ψ = Ext S 0 (e ; ϕ) from S 0 of the scalar field ϕ in the direction e as An equivalent variational characterization of this extension operator is We can define a norm on the surface-based data using this extension, namely, ψ E 0 . Typically, considering the following linear observation operator we can use this norm in the observation space. Accordingly, observer terms in the form H * (z − H u) will give in a variational setting where the second expression is obtained by using the characterization (9) when observing that on S 0 Therefore, we have in this case Now, when dealing more generally with a discrepancy operator J( u, z) pertaining to the displacements on a surface S 0 , we will generalize this strategy by using the observer correction given by obtained by directly substituting in (10) J( u, z) for (z − H u), and −D u J( u, z)-associated with a vector field on S 0 -for e since in the linear case we have H u = u · e on S 0 . Of course, this also easily generalizes to a measurement made of a collection of such surfacebased data associated with N S surfaces (S i 0 ) N S i=1 and associated vector fields e i , for which the correction will be given (in the linear case) by

Tag planes
We consider data consisting of a set of N P tag planes T = N P i P i deforming over time. Following the original ideas of [6] the discrepancy between the model and the data will be measured using the signed distances between the tag planes and the corresponding model data. These model data are deforming surfaces obtained by applying the model displacements to the initial configuration of the tag planes. Let us then denote by T 0 = N P i P i 0 the set of tag planes in the reference configuration, mapped by the estimated trajectory u to T = N P i P i . For any point in a model tag plane x = ξ + u(ξ) ∈ P i for some ξ ∈ P i 0 , we can compute the signed distance to the corresponding target tag plane where x is a point on the ith model tag plane P i , P i x is the projection of this point on the corresponding observed tag plane, and n P i is the normal of the observed tag plane at the projection point-see Fig. 2. The discrepancy operator is then the application mapping the displacement field to this collection of (scalar) distance fields defined over the planes of T 0 . When differentiating with respect to the displacement field we have  Example of a tag grid and apparent displacements extracted from real tagged-MR images. Apparent displacements obtained from inTag plugin of OsiriX software [48] Hence, the application of the above-described strategy gives an observer that follows the mechanical system of equations (3), except for the first equation modified intȯ

Tag grids
We now consider the data in the form of a collection of tag lines deforming within the set of (2D) image slices-see Fig. 3a for an example. We thus assume that we have N P such lines L ij N P i=1 in each 2D image I j , with 1 ≤ j ≤ N I . We cannot directly design the discrepancy operator based on the corresponding model lines, since displacement fields are not well-defined along lines in the variational space. To circumvent this difficulty, we again consider the tag planes in the model and project each point of the planes onto the neighboring image slices. Denoting by I j the Euclidean orthogonal projection onto the image I j , we compute the signed distance of the projected point to the corresponding tag line within each image-see Fig. 4 Then, we can interpolate the signed distances thus-obtained in the various images concerned-which provides interpolated distance fields over the model tag planes as a discrepancy operator, namely, for each plane P i 0 , where J (j) denotes the interpolation operator. When differentiating this expression, we have but we also have a contribution coming from the interpolation operator derivative. Since this interpolation only depends on the coordinate of the point considered along the axis orthogonal to all image slices, denoting by J (j) the derivative with respect to this coordinate, a straightforward computation finally yields where n I denotes the direction vector of the orthogonal axis. Note that when considering e.g. linear interpolation, the derivative J (j) is directly given by the finite difference expression computed between the two adjacent planes. This finally gives for the observer correction equatioṅ

Filtering 2D apparent displacements
Finally, we consider the measurements corresponding to apparent displacements obtained from the processing of 2D tagged images-see Fig. 3b for an example. In the literature we can distinguish two main corresponding extraction methods. A first manner consists in tracking the tag pattern directly in the image plane. For instance [49,50] consider an optical flow methodology that takes into account the fading of the tag pattern during the acquisition process. Another solution proposed by [41] is to perform non-rigid image registration. A second family of methods consists in working in the frequency domain. The most popular method is the HARP technique-see [51,52]-which tracks the phase of the tag pattern. Following this trend, recent works proposed by [42,53] use the Gabor filter to obtain a better estimation of local deformations in late systole-which appears to be a slight limitation of the HARP methodology.
In any case, we can assume that this pre-processing step enables us to track-throughout the dynamic sequence-the intersections of material fibers originally orthogonal to the image planes. This holds e.g. for such fibers corresponding to the intersection of tag planes, but 2D-tag processing techniques generally provide a (2D) displacement field all over the image slices-as depicted in Fig. 3b. Note that the material point located at the intersection between the fiber and the image changes over time due to through-theplane motion. Hence, the measurement is not a material displacement, see Fig. 5. This induces serious complications in the exact form of the tangent observation operator DH, therefore we will use an approximate form based on a small displacements assumption. With this assumption, the observation operator reduces to the components of the material displacements tangential to the image plane. In this case, the measurement is a two- where e denotes the projection onto the plane orthogonal to e, plane in which ϕ is assumed to lie. Finally, the correction equation for the observer readṡ where n I denotes the normal to the image planes, and the innovation term z − H( u) will be computed-see "Generating apparent displacements" section-based on the actual tracking of material fibers, i.e. without small displacements assumption.

Time discretization of the state observer
In this section we address the issue of the time discretization of the observer. Here, we focus our effort on ensuring that the dissipative behavior of the (time discrete) estimation error dynamics is preserved, up to some consistency terms inherent to any discretization. A specific numerical time scheme based on a mid-point scheme was proposed by [6] for similar observers. Our present approach, however, differs in the sense that the time-discrete observer is built on a prediction-correction paradigm. Consequently, the prediction part-in practice, iterations of the direct model-and the correction parti.e. the action of exploiting the discrepancies between the data and the model-can be managed in separate ways.
In order to incorporate the observation operators described in the previous section, we consider the case of nonlinear operators. More precisely, neglecting the observation noise, we assume that the observations are obtained by where H(·) is a nonlinear and sufficiently smooth observation operator. Following [6] we propose to define the time-discrete observer as where we denote by y n+1 − and y n+1 + the prediction and correction steps respectively. In the correction relation (19b), we use y e + an extrapolated trajectory, and DH e + the tangent operator of the observation operator evaluated at y e + , i.e. DH e + = DH( y e + ). While other choices are possible-see [6]-the most simple approach is to set y e + = y n+1 − . Even though fewer theoretical results can be obtained in the case of nonlinear observation operators, this numerical procedure is based on a linearization argument. This leads, after a local analysis around the trajectory used in the linearization, to a dissipative behavior similar to that of a linear problem. To fully understand this crucial aspect, we derive in Proposition B.2 of Appendix B the energy estimate satisfied by the estimation error, in the fully linear case. Building on this first result, we deduce in Proposition B.4 of Appendix B a similar relation satisfied by the linearized estimation error. Due to the linearization procedure proposed in (19b)-compared to an explicit scheme-this estimate shows desirable dissipation properties, up to two (natural) source terms: a first consistency term emanating from the time discretization process, and a second term due to the linearization process. Hence, assuming that the initial condition is reasonably close to the target initial condition, the latter will have little influence on the overall stability of the numerical scheme.
In the case of discrepancy measures-see "Filtering tagged-MR planes and grid" section for practical examples-the observations and the real trajectory are linked through the implicit relation J y(n t), z n = 0.
In this case, the time-discrete observer can be directly inferred from system (19a)-(19b) y 0 + = y 0 The corresponding linearized estimation error satisfies the estimate in Proposition B.4 of Appendix B, but the operator DJ e + = DJ( y e + , z n+1 ) appears instead of the tangent of the observation operator.

Characterization and iterative resolution of the correction step
We now introduce a fully discrete model-solved in practice-by considering the vectors of degrees of freedom associated with a standard finite element spatial discretization. We denote by capital letters the vectors of degrees of freedom, and by italic operators the matrices associated with the functional operators used until now. For example, we denote by U ∈ R N dof the vector of degrees of freedom-of dimension N dof -associated with the function u h defined in the finite-element space Y u h ⊂ Y u . We also define a state vector as the concatenation of the displacement degrees of freedom and the velocity degrees of freedom Y = U V . This final space discretization procedure enables us to define a fullydiscrete transition operator F n+1|n , such that the prediction step (19a)-or (20a)-can be rewritten as where F n+1|n represents the time-discrete flow from time t n to time t n+1 . Concerning the correction step, some specific elements need to be addressed due to the presence of the adjoint of the tangent observation operator in (19b)-or of the discrepancy operator in (20b). This aspect is particularly important-as extensively detailed in [6]-when dealing with displacement-based observations. To that end, we define for all (u h, and the mass and stiffness matrix respectively, and we consider the associated norm In the same manner, we denote the linear operator after space discretization by H, and by DH the tangent of a nonlinear observation operator. By a slight abuse of notation we use the same notation H-or DH-when applied to U or to Y , despite the fact that it corresponds in this latter case to H 0 -or DH 0 . Concerning the observation norm, we consider the matrix S computed through the extension operators as detailed in the previous section. For example, when considering 3D displacements, let us consider Z h ⊂ Z a suitable discretization of the observation space. The (linear) observation operator boils down to an interpolator between the mesh and the observation subdomain ω 0 , whereas S is defined for all (ψ h , 1 We recall that, in this particular case, the extension operator is defined by (4). This extends directly to tag planes (or tag grids) and to apparent displacements using the definition of the extension operators (8) and (16) respectively. Taking into account this spatial discretization, we can write the correction step (19b) after space discretization as where DH e + is the tangent of the observation operator around the extrapolated trajectory Y e + . To simplify the presentation, we set Y e + = Y n+1 − , so that the correction step can be expressed in a more practical form as where S t = tS. We remark that, in a fully linear setting, (22) corresponds exactly to the Best Linear Unbiased Estimator (BLUE) associated with the observation Z n+1 , the observation covariance S −1 t , the a priori state Y n+1 − and the a priori state covariance γ N −1 [15]. Alternatively, it can be interpreted as solving the following minimization problem where · N and · S t are the norms induced by the matrices N and S t respectively. From this alternate characterization, we see that the correction step corresponds essentially to a linearized version of a least-squares minimization around the a priori state Y n+1 − . Note that similar characterizations can be deduced from (20b) in the case of discrepancy operators. The complete system (21)-(22) can be seen as a prediction-correction discrete-time sequential estimator as is the case for the discrete-time Kalman filter [15,54], but here the a priori state covariance remains constant equal to γ N −1 . Therefore, whereas the discrete-time Kalman filter is not computable for systems arising from PDEs, our filter is, since N is sparse.
Let us now give some additional methodological key points to solve the correction step (22) in a very effective way. We remark that which proves that N −1 is a good preconditioner to solve Eq. (22) with an iterative solver. This reveals to be very helpful as typically we do not want to store the operator DH e + S t DH e + . Indeed, with an iterative solver, we only need to be able to compute for any vector Y quantities like N −1 Y and (N + γ DH e + S t DH e + )Y . Using an iterative solver is particularly effective in the case of apparent displacements where the observation space is the concatenation of the set of image planes having a potentially high resolution, yielding a very dense operator S.
Finally, a practical difficulty may arise in our specific case where quantities of the form DH e + S t DH e + Y may be cumbersome to compute because of the choice of the observation norm S, obtained from various extension operators. Authors in [6] demonstrated that a judicious approximation of the extension can be found by considering a penalized minimization problem. For instance, in the case of 3D displacement measurements, the computation of the adjoint ψ = H * ϕ = Ext ω 0 (ϕ) for ϕ ∈ Z can be approximated, after space discretization, by solving min ∈R N dof where ε > 0 is the (small) penalization parameter, and are vectors of degrees of freedom, H is an interpolation operator between the mesh and the subdomain ω 0 , and M obs is the matrix associated with the L 2 -norm on the observation space From the penalized form (23) of the extension operator, we derive a more convenient expression of the observation norm S. To do so, we start by remarking that, for all where we have used the orthogonality property (5) to obtain the last relation. Then, replacing the extension operator by its penalized approximation (23) yields from which we retrieve the relation given in [6], Similarly, when dealing with tag planes, tag grid or apparent displacement, one can approximate the extension operators (8) and (16) using a penalization strategy.

Time discretization of the joint state and parameter observer
Having defined the discrepancy operator and designed the state observer, we can now consider the additional stage of parameter identification through the state and parameter observer introduced at the end of "Principles of sequential estimator design" section.
We should now define the adequate discretization of System (2) compatible with the discretization of the state observer already defined in (20a)-(20b). To that purpose, two discretizations are available in the literature. The first one is based on the fact that (2) corresponds to applying a continuous-time Reduced-Order Extended Kalman Filter (RoEKF) to the parametric space, hence a proper discretization is clearly the prediction-correction scheme defined by the discrete-time RoEKF [18]. The second one proposed by [55] is not directly an exact discretization but rather an extension at the discrete-time level. In fact, the parameter dependency makes the joint state and parameter system nonlinear even if the state dynamics is linear. Therefore, the RoEKF filter on the parameters is only an approximate optimal filter. Other choices of approximate reduced-order optimal filter can therefore be used when available. It is typically the case of the Reduced-Order Unscented Kalman Filter (RoUKF) derived by [55]. This filter replaces at the discrete-time level the tangent computations in (2) by finite difference computations which appear to be better adapted to large nonlinearities. In addition, there is no need to implement the tangent operators, as instead the original dynamics and observation operators are applied on so-called "sampling points". This algorithm thus combines accuracy and computational efficiency, and it has already been successfully applied in real biomechanical applications, indeed, see [7,21,56,57]. Moreover, it reduces to the Reduced-Order Kalman Filter upon linearization, which allows to validate its stability with an error linearization study as achieved by [55]. Indeed, both algorithms reduce to the reduced-order Kalman filter after linearization. For completeness we here recall the complete algorithm in our case, before proceeding to the results section. Following [55,58], for θ discretized in R r , we introduce r + 1 so-called unitary simplex sampling points I [i] in the space R r and the associated weights α i with the following rules i.e. with zero mean and unit covariance so that, at each time step, the sampling points can be generated around the estimated values based on the actual covariance estimation. Given an adequate sampling rule, we store the corresponding weights in the diagonal matrix M α and precompute these unitary sigma-points. We also denote by M t = tM obs , and by [I [ * ] ] the matrix concatenating the (I [i] ) 1≤i≤r+1 vectors side by side, and similarly for other matrices aggregating some particle vectors. We then perform at each time step 4. Parametric correction:
Using the calibration strategy of [25]-based on a reduced modeling approach-we propose a complete set of model parameters, see Table 1, for which direct simulations are in good agreement with standard values of physiological and mechanical indicators-see Fig. 7. For the parameters describing the infarct we choose (θ K , θ) = (1, −1) in the septum, and (θ K , θ) = (0, 0) otherwise.

Synthetic data generation
Generating the synthetic measurements required to assess our proposed tagged-MRI data assimilation strategy is an intricate procedure in itself. In this section we discuss the various methodological steps used to build, from the biomechanical heart model, the tag planes, the tag grids and the 2D apparent displacements.

Generating tag planes and tag grids
A natural idea to build the set of tag planes from a direct simulation is to produce, in the deformed configuration at marking time, a set of two-dimensional triangular meshes associated with the planes and to consider the nodal interpolation operator between the model tetrahedral mesh and the set of tag meshes. Denoting by J P i m the interpolation operator of a single tag plane, its displacement u P i m is given in this context by where u m is the model displacement at marking time t m , ξ m denotes the associated deformed coordinate, namely, ξ m = ϕ(ξ, t m ), and m = ϕ( 0 , t m ). This leads to significantly irregular displacements near the intersections between the model boundary and the tag planes. One way to circumvent this limitation is to consider the tag planes as made of an elastic material and to regularize the interpolated displacement using an appropriate elastic model. However, as the geometry at hand is two-dimensional, a shell model would be required. To simplify this task-which is, in essence, an issue of data regularizationwe consider a set of elastic (3D) tag layers In practice, each tag layer is built so that P i m ⊂ V i m . Hence, the displacement of a tag plane P i m is derived from where u V i m is the displacement of the tag layer V i m , verifying The procedure described in (28) is in fact the extension-in the sense of "Filtering tagged-MR planes and grid" section-of the interpolated displacement, namely, In (28) we denoted by σ V the Cauchy stress tensor describing the tag layer material. In practice, a relevant choice is a linearization around a given trajectory of the heart material since, within the image, the tag pattern follows the heart material points. We show in Fig. 10a an example of synthetic tag planes extracted from a direct simulation. As far as the tag grids are concerned, they are obtained by clipping the tag plane meshes with the image planes. Figure 8 illustrates the complete procedure of construction of a tag plane and of several tag lines, while Fig. 10b gives an example of extracted synthetic tag grids.

Generating apparent displacements
Once the tag grids are created, the apparent displacement field can be approximated by tracking the displacements of the tag lines intersection points. More precisely, at marking time, we compute the intersection point of every tag lines ((red) crosses in Fig. 9a). During the simulation, as the tag lines deform, we track the displacements of the intersection points, leading to the (green) vectors in Fig. 9b. Once the displacements of the intersection points are computed, a global apparent displacement field of the image plane is produced ((blue) vectors in Fig. 9c) by standard interpolation. We show in Fig. 10c an example of synthetic apparent displacements extracted from a direct simulation.

Discrepancy measure in practice
The tagging process is not performed in the reference configuration-never observed in reality-but at the previously mentioned marking time. For this reason, the tagging pattern is necessarily built over an already deformed configuration. Hence, the information obtained from a set of tagged-MR images should be considered of Eulerian nature. However, assuming the displacement at marking time u(ξ , t m ) is given, we can circumvent this difficulty by introducing this additional information in the filtering procedure. From an algorithmic standpoint, the discrepancy measure is computed as follows: 1.
[offline] Build at marking time the set of tag planes from the deformed configuration (by u(ξ , t m )) to the tag planes [online] From the estimated displacement, deform the tag planes [online] From the estimated (deformed) tag planes N P i P i compute the innovation terms appearing in (12) (planes data), in (14) (tag grids) or in (17) (apparent displacement).
In our context of synthetic data assimilation, we directly provided the displacement u(ξ , t m ). In real cases, the task of estimating the displacement at marking time could be carried out using, for instance, the segmentation of the endo-and epicardium of the left ventricle-obtained typically from cine-MR images. In any case it requires another source of information on the system and this points out a certain limitation of the tagged-MRI data for estimation purposes.

Spectral analysis of the observer associated with (3D) tag planes
In this section we discuss the observer built using a set of tag planes that can be decomposed into three distinct families. Each tag plane in a given family shares-at marking time-the same orthogonal direction and, only for the sake of simplicity, we assume that the three directions are orthogonal. This type of data set may be referred to, in what follows, as (3D) tag planes. As already mentioned and emphasized in Appendix B, the quality of the state filtering procedure can be assessed by the amount of damping we introduce in the otherwise conservative or weakly damped system. For this reason we propose to conduct a numerical study of the spectra of the operators driving the target dynamical system and the estimation error dynamical systems. For simplicity, we consider a linear elastic model, typically obtained from the linearization around the null trajectory of the calibrated passive cardiac model. To facilitate this analysis we also decrease the natural viscosity of the target system. Hence, we consider the solutions of the following spectral problem which corresponds to the operator without filter and where, additionally to the stiffness K and mass M matrices, we denote by C the damping matrix obtained after spatial discretization. For the case of complete displacement observer we consider the spectral problem where H corresponds to the identity matrix or an interpolation operator. Note that the form of the stabilization operator is the one derived in (24) where the extension operator has been approximated by its penalized counterpart. For the (3D) tag plane observer, the where the discrepancy operator corresponds to the one detailed in (13). Using the optimal criterion on the gain γ provided by [6], in Fig. 11 we show the spectra obtained for the operators without filter, with complete displacement feedback and with (3D) tag planes. In these plots we also vary the spacing between two consecutive tag planes from 8 mm, to 3.5 mm and 0.25 mm. In the three situations we observe that for low frequencies the observer using tag planes acts as the direct displacement observer. For coarse tag patterns we naturally observe that the higher frequencies are less stabilized than with the direct displacement observer. This phenomenon disappears as the tag pattern becomes denser, and the overall efficiency of the tag planes observer tends towards the full displacement observer. This result is of major importance since it comforts the intuition-and hopes-that tagged-MRI could indeed provide information on intra-myocardial deformations. Even though a more thorough mathematical analysis should be carried out to corroborate these results, a first explanation for these outstanding performances is that the tangent of the discrepancy operator is the concatenation of the normal vector fields of the tag planes, which in the case of 3D tag planes span all directions in space.

Estimation results
Since we have evaluated-through spectral analysis-the state estimation capabilities of our proposed observer, we can now assess the joint state and parameter estimator by identification of the infarct location and intensity presented in the beginning of the results section. We distinguish two situations. First we provide the exact location of the infarcted region but not the values of the parameters characterizing the pathology. Hence, we have two separate regions, for each of which we identify a contractility parameter and the main passive stiffness parameter. This case may be referred to as the 2-Regions case. Secondly, we use the AHA-Regions-roughly corresponding to the territories of the coronary arteries, see [59]-to partition the left ventricle and retrieve both active contractility and passive stiffness parameters in each AHA-Region. In this case, the infarct location will be inferred from the parameters spatial variation over the AHA subdivision. In all cases, for the sake of simplicity, the reference configuration and the intra-cavity pressures are assumed to be known and the initial condition is defined by solving an equilibrium state with the lowest pressure sustained before the atrial contraction. We point out that, since the stiffness is globally modified between the target system and the observer, an error in the initial condition will be introduced during the estimation.

The 2-Regions case
We start with the simpler 2-Regions case and we only seek to assess the observability provided by the data. We choose a fine spatial distribution of the tag planes by setting the space between two consecutive tag planes to 3.5 mm. Denoting by M T the surface mass matrix computed on the set of tags, we define the measurement observation norm S by setting, in (23) The parameter m represents the square of the standard deviation of the discrepancy measure. In the perspective of only assessing the method capabilities, the observations are extracted from the direct simulations-as explained in "Synthetic data generation" section-with a high temporal resolution of 1 output every 25 simulation time stepsset in our simulations to 2.5 · 10 −4 s-and no noise is added. Therefore, following [7] we rescale m = t obs t m obs = 25 m obs , and set a high confidence in the observations with m obs = (0.65 mm) 2 . The results are presented in Figs. 12, 13, 14 and 15, where the dashed curves visualize the estimated standard deviation. In Fig. 12, we consider the most optimistic configuration where threedimensional tag planes are available, whereas in Fig. 13 we consider a more realistic configuration where only two directions of tag planes-here a short axis grid only-are available. Then, in Fig. 14 we proceed with the corresponding bi-dimensional grid as described in "Synthetic data generation" section. We thus rely on the grid-based observer discrepancy. Finally,we present in Fig. 15 the results where extracted 2D apparent displacements are defined as the available measurements. In this particular case, since the innovation corresponds to the comparison of two vector fields defined within the image planes, we take the observation norm as a piecewise constant mass matrix in the image domain. The behavior of the estimation procedure is very similar for all types of processed data considered, and the parameter values are accurately estimated in the two regions, both for the active (contractility) and passive (stiffness) parameters. It should be noted, in particular, that the estimations produced based on 2D tagged data-namely, with tag planes and grids, and apparent displacements-are as effective and accurate as those obtained with 3D tags.

The AHA-Regions case
We now consider a configuration that can be encountered in a more realistic context. First, we only generate from the direct simulation approximately 20 "synthetic" images in the cardiac cycle from which we extract 2D tag planes. Then, for the estimation, we set up an AHA-based partition, which defines 17 regions, namely, the 16 AHA regions and the remaining part of the heart. The evolution of the joint state and parameter estimation is presented in Figs. 16, 17 and 18. The state convergence is demonstrated through the evolution of the volume curves and P-V loop plots, whereas the evolution of each parameter is presented in Figs. 17 and 18. The results are divided into three groups of parameters associated with the three different long axis elevations of the AHA partitions, namely basal (regions 1-6), intermediate (regions [7][8][9][10][11][12] and apical (regions [14][15][16]. We visualize in Fig. 19 the final parameter identification diagram in a bull's eye representation as used in a clinical context. In addition, we recall the identification that would be obtained when only exploiting cine-MRI segmentations of the endocardium and epicardium as presented in the previous work of [7]. This directly illustrates the identification benefits obtained with the tagged-MRI measurements, in particular concerning the passive stiffness parameters.

Discussion
First, concerning state estimation, our results confirm the remarkable effectiveness of the strategy previously introduced by [6], and the excellent adequacy of both tagged-MR as an imaging modality of choice for estimation purposes, and of our herein-proposed strategies for incorporating such data via specifically designed discrepancy operators. Indeed, the above spectral analysis gives a very clear indication as to how fast state estimation errors are being damped when using this estimation chain. In particular, we have observed the convergence of the spectrum-with respect to tag spacing-towards that of the observer with full 3D observation, while the spectrum obtained with coarse tags shows that standard tag spacing is amply sufficient to obtain uniform damping rates. This is also confirmed with the results obtained in the joint state-parameter estimation trials, in which mechanical indicators are effectively and accurately retrieved, recall Fig. 16.
As regards parameter estimation, the additional estimation stage provided by RoUKF filtering-combined with the state observer-also shows very good performance. In the 2-Regions estimation setup, in particular, both active and passive parameters are very accurately estimated, and in a very short time as soon as the parameters concerned become observable in the type of behavior that is encountered along the cardiac cycle. Namely, passive parameters are mostly observable during the-rather short-initial diastolic phase associated with atrial contraction, while of course contractility parameters can only be revealed once the electrical activation actually starts. Note that this 2-Regions setup gives a realistic strategy in clinical perspectives, as cardiac MR performed for infarct diagnosis frequently includes late-enhancement sequences, which can be segmented to provide the desired subdivision into healthy and diseased regions. In case late-enhancement images (or the associated segmentations) are not available, or when additional concurrent localization information is desired, estimation can be performed based on the AHA subdivision. The corresponding estimation results exhibit the same general features as with two regions-namely, rather fast convergence during diastole and systole for passive and active parameters, respectively -albeit as expected the estimation is less accurate for each individual parameter. Nevertheless, active parameters are still quantitatively retrieved, and passive parameters are quite discriminately detected within the infarcted region, and much more so than with estimation based on Cine-MR. Of course, fundamental identifiability issues are of concern in this multiple parameter estimation context, and we can expect that identifiability would be improved-hence estimation would be more accurate-when using segmented Cine sequences in addition to tagged images in the estimation procedure.

Towards clinical applications
Despite its efficiency with synthetic data, our data assimilation strategy still needs to be validated using real clinical data. To address this challenge, it remains to complete the proposed approach with the following requirements: 1. Obtaining patient specific geometries and fiber orientations. Even if segmenting the ventricles is a well studied problem, it remains an issue [3,60]. It is even more the case for obtaining fiber orientation maps, directly from images [61], or by combining geometrical models [62] whose parameters could also be assimilated [63]. Moreover, with large deformation systems such as the heart complete bio-mechanical model, we ultimately need a reference configuration which needs to be inferred from the images at hand, typically using inverse methods [64]. 2. Registering the tagged-MR sequence with respect to the initial geometry, hence to the reference configuration by transitivity. 3. Including blood pressure measurements. Note that we should either rely on direct but invasive measurements in the cavities or on estimating these blood pressures from non-invasive distal pressure using typically signal processing [65,66] or data assimilation.
4. Specifying the boundary conditions, or using additional measurements to estimate them. In our modeling choice, viscoelastic boundary conditions are defined on n which needs to be geometrically designated. Moreover, this simplified model should be ultimately jointly estimated with similar data assimilation methods [21].
One tremendous difficulty is that each error in the steps listed above can introduce irreducible errors in the estimation procedure. This was typically our experience in [7,21,67] where we discussed how boundary conditions modeling errors can pollute the resulting estimation in a "clinical" context. Therefore, even if our method paves the way for carefully integrating tagged-MRI measurements in a data assimilation strategy, the complete model-data fusion pipeline in clinical cardiac applications necessitates to combine multiple technologies and expertise from modeling to data processing and associated registrations [3,5].

Conclusions
We have proposed specific methods for integrating tagged-MR sequences in a data assimilation framework with a beating heart model. Tagged-MRI represents the "gold standard" in cardiac imaging, and great benefits are expected from using the corresponding rich kinematical information for performing the joint state-parameter estimation of the system, and of various modeling parameters of high potential value in terms of clinical diagnosis assistance.
In this data assimilation framework, a crucial ingredient lies in the adequate formulation of a discrepancy operator to compare the model and the data. We have considered several options, based on: (1) extracted 3D displacements; (2) tag planes in the 3D volume; (3) tag grids in 2D slices; (4) apparent displacements in 2D slices. In practice, the specific choice of discrepancy operator could be based on the type of tagged sequence and on the available corresponding post-processing tools, albeit our unified framework also allows detailed comparative assessments. For the purpose of state estimation, each definition of discrepancy operator was accompanied by the formulation of an adapted filtering operator.
We have also proposed well-adapted discretization strategies. As regards time discretization, in particular, a two-step "prediction-correction" type algorithm was designed for the proposed estimation systems, allowing to completely dissociate the operations related to the model and those performed for estimation purposes, e.g. with two differentcoupled-software codes. This is very valuable from a software architecture perspective, and in particular makes the estimation strategy compatible with modular concepts such as those underlying the Verdandi data assimilation library [68].
Mathematical analyses have been provided at the various stages of construction, to substantiate the convergence of the overall observer strategy based on a simplified-albeit illuminating-linear example, and also to assess the effects of discretization procedures.
Finally, some detailed numerical assessments of the overall estimation framework have been performed, based on synthetic data produced by a reference cardiac simulation representing the behavior of an infarcted heart in a realistic manner. The assessment results show that state estimation is extremely effective, while the performance of parameter estimation depends on the specific estimation objectives, as can be expected from the point of view of observability. In particular, when the diseased region is pre-determined prior to estimation, active and passive parameters are very accurately and quickly retrieved in the infarcted and healthy regions. When the more challenging objective of estimation in an AHA subdivision is considered without any prior on the diseased region, the convergence of each individual parameter value is less accurate, but the overall distribution of parameters is very adequately retrieved, allowing for effective localization and quantitative assessment of the disease. This provides a great improvement over similar estimation based on using Cine-MR alone, which only gives adequate results for active parameters.
All major ingredients are thus in place for using this methodological framework in a patient-specific context with actual data, which is of course a most natural perspective of this work. Addressing this challenge will imply combining different measurements and image modalities in order to estimate patient-specific modeling components. Other perspectives concern the consideration of alternative discrepancy operators, such as with the formalism of currents [69][70][71] which would allow to dispense with using sophisticated image post-processing tools on tagged-MR as a prerequisite for data assimilation.
In a semi-group theory context we introduce the semi-group generator A ∈ L(D(A), Y) with and we can prove that (29) admits a classical solution Moreover, the operator A is skew-adjoint, implying that corresponding to the energy balance on the system (29) with E = 1 2 y 2 Y . Using the semi-group theory, we rewrite the dynamics of the model in the abstract state-space forṁ y = Ay + R, y(0) = y 0 + ζ y .
Concerning this model, we assume that we have at our disposal some measurements of the displacements. We introduce the observation operator H = (H |Y u 0) such that where ω 0 is an internal open subdomain of 0 and Z = H 1 (ω 0 ) 3 . Here H and H |Y u as H = (H |Y u 0) correspond to the same operators with different input spaces. Using the extension ψ = Ext ω 0 (ϕ) defined by we can prove the following property.
Proof The proof is a simple extension of the property proven by [73] for scalar equations. Let ϕ be an element of Z. The only difficulty lies in proving the norm equivalence with ϕ 2 H 1 (ω 0 ) 3 . First, we have with C p given by the Poincaré inequality and C k given by Korn inequality and a bound C a on the elasticity tensor. Conversely, by continuity of the extension on 0 \ω 0 with respect to the data, there exists a constant C d > 0 such that for any ψ = Ext ω 0 (ϕ) we have Hence, denoting by C t the constant arising from the continuity of the trace operator, we have which completes the proof.
It is now possible to define the adjoint of the observation operator.
Proposition A. 2 The operator H is bounded from Y to Z and H * is given by Proof Let us first prove that H is bounded. We consider ψ ∈ Y u and ϕ such that ϕ = H |Y u ψ. We have directly, from norm equivalences, By the variational characterization of the extension (5) we have 0 ε(Ext ω 0 (ϕ)) : A : ε(Ext ω 0 (v |ω 0 )) d = 0 ε(Ext ω 0 (ϕ)) : A : ε(v ) d , and H * is given by We can now define the observer by the dynamics which in strong form reads which converges to the solution of (29) under the observability condition given by the next theorem, see e.g. [74] for a proof.
We can now make explicit the specific observability condition in our configuration that will allow us to invoke Theorem A.3 of Appendix A.

Theorem A.4 If there exists a constant C st and a time T such that every solution of
then, in the absence of observation error, the observer given by the dynamics (31) converges to the solution y ref of (29) such that z = Hy ref .
Proof which, from Theorem A.3 of Appendix A, converges exponentially to 0 for every initial condition when the observability condition (33) is verified.
Following [75], we define the elastic geometric control condition: The elastic geometric control condition is satisfied if every combination of pressure (P) and shear (S) waves ray encountersin the sense of [75]-the subdomain of observation.
Readers may refer to [75] for a complete description of such rays. This condition generalizes to the vectorial case the so-called geometric control condition (GCC) introduced by [76], allowing to control any solution of the acoustic wave equation from the observations of the time derivative of the wave in a subdomain.
Theorem A. 5 The observability condition of Theorem A.4 holds on a subdomain ω 0 as soon as the elastic geometric control condition is satisfied withω 0 and dist( 0 \ω 0 ,ω 0 ) > 0.
Proof For technical reasons, we assume that the elastic geometric control condition is satisfied for an observation domainω 0 slightly smaller than ω 0 , namely, withω 0 ⊂ ω 0 and dist( 0 \ω 0 ,ω 0 ) > 0. We first recall the classical observability result when the velocity is observed. In fact there exists a constant C st and a time T such that every solution of (29) satisfies the observability condition withT = T − δ for δ > 0 sufficiently small, as soon as the elasticity geometric control condition is verified in the time interval [0, T [ [75]. Following what was already done for acoustic waves by [73] we will use a property of equirepartition (over time) of the total energy localized within the observation subdomain between the kinetic and elastic contributions to infer (33) from (34).
Let ψ ∈ C ∞ c ( 0 ) be a cutoff function satisfying

Linear model and linear observation operator
To start with, we assume that the dynamical system satisfied by the target trajectory is given by ⎧ ⎨ ⎩ẏ (t) = (A + ηV)y(t) y(0) = y 0 + ζ y , where A is a skew-adjoint operator, V is a self-adjoint and semi-negative operator and η ≥ 0 is a viscosity coefficient. Note that (37) can be interpreted as a linearization of (3) around the stress-free configuration. Additionally, we consider a linear observation operator and we neglect, for simplicity, the observation noise. Denoting by t the (constant) time step of the numerical procedure, the time-discrete observations read z n = Hy(n t). The prediction-correction scheme for the observer reads ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Relation (38a) corresponds to the prediction step, with the operators driving the target system, while (38b) is the correction step. Defining the discrete estimation error from the correction step y n + = y(n t) − y n + , and associating a prediction error with y n+1 − = y((n + 1) t) − y n+1 − , we can determine the time-discrete dynamics satisfied by the estimation error. Proof (41b) is directly inferred from the definition of the prediction estimation error and using (38b), namely, y n+1 − = y((n + 1) t) − y n+1 − = y((n + 1) t) − y n+1 + + γ tH * H y n+1 We now have to work our way to (41a). First, we remark that y n+1 − − y n + t = y((n + 1) t) − y(n t) t − y n+1 − − y n + t .
In (45) we see the effect of the correction step (38b) leading to some dissipation terms brought by the observation operator. The expression is the abstract and discrete version of expression (7), perturbed with natural consistency terms.

Linear model and nonlinear observation operator
As a second step, we consider the case of a nonlinear observation operator, so that the observations are obtained through (18), while the model operator remains linear. In this case, the related prediction-correction time-scheme is This entails (48b) since from equation (47b) and the linearization of H(y((n + 1) t)) around the extrapolated trajectory y e + we have y n+1 − = y((n + 1) t) − y n+1 + + γ tDH e + * z n+1 − H( y e + ) − DH e + ( y n+1 + − y e + ) = (1 + γ tDH e + * DH e + ) y n+1 + + tλ n .
All the other computations previously presented to prove Proposition B.1 still hold, so that we obtain the dynamical system (48a)-(48b) satisfied by the estimation error.
We can now establish the energy estimate associated with (48a)-(48b).

Proposition B.4
The norm of the estimation error in the case of a nonlinear observation operator satisfies the following estimate Proof We remark that the system satisfied by the linearized estimation error (48a)-(48b) can be obtained from system (41a)-(41b) by replacing the observation operator by its tangent around the linearization trajectory, and by adding the linearization term λ n in the correction step. Hence, by following the demonstration of Proposition B.2 of Appendix B, one can obtain estimate (49).
Finally, let us remark that the extension of these results to the case of a nonlinear model does not entail specific challenges, within the context of the presented analysis. After linearization, one can expect to produce another linearization source term, this time in the prediction step (48a)-instead of the correction step. Hence, the dissipation property of the time-discrete observer still holds, up to this additional linearization term.