Research article
Open access
Published: 09 January 2021

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

Advanced Modeling and Simulation in Engineering Sciences volume 8, Article number: 2 (2021) Cite this article

2893 Accesses
4 Citations
Metrics details

Abstract

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 internally—in 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.

Position of the problem

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 ${{\mathrm {y}}}^{{\text {ref}}}(t),\ t \in [0,+\infty )$, belonging to the so-called state space ${\mathcal {Y}}$ and solution, in our case, of a—possibly infinite-dimensional—dynamical system summarized in the state space form

$$\begin{aligned} {\dot{{{\mathrm {y}}}}}^{{\text {ref}}} = {{\mathrm {F}}}({{\mathrm {y}}}^{{\text {ref}}},t), \end{aligned}$$

with an uncertain initial condition

$$\begin{aligned} {{\mathrm {y}}}^{{\text {ref}}}(0) = {{\mathrm {y}}}_0 + \zeta _{{\mathrm {y}}}, \end{aligned}$$

where ${{\mathrm {y}}}_0$ is a known a priori and $\zeta _{{\mathrm {y}}}$ is the uncertain part in the initial condition. Therefore, any simulation of ${{\mathrm {y}}}$—based on the discretization of the dynamical system—starting only from ${{\mathrm {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 ${{\mathrm {z}}}$ these measurements—also called observations and belonging to the observation space ${\mathcal {Z}}$—which are assumed to be generated by a mapping ${{\mathrm {H}}}$ on the real trajectory, up to additional measurement errors

$$\begin{aligned} {{\mathrm {z}}}= {{\mathrm {H}}}({{\mathrm {y}}}^{{\text {ref}}},t) + \chi . \end{aligned}$$

The observer denoted by ${\widehat{{{\mathrm {y}}}}}$ is a system that starts from the only part known in the initial condition—namely ${{\mathrm {y}}}_0$—and uses in time the available measurements ${{\mathrm {z}}}$ to generate a trajectory ${\widehat{{{\mathrm {y}}}}}(t),\ t \in [0,+\infty )$ that converges to ${{\mathrm {y}}}^{{\text {ref}}}$ as fast as possible. Therefore, simulating ${\widehat{{{\mathrm {y}}}}}$ instead of ${{\mathrm {y}}}$ from ${{\mathrm {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

where ${{\mathrm {G}}}$ is called the gain operator, also referred to as filter. The dynamics of ${\widehat{{{\mathrm {y}}}}}$ is corrected when a discrepancy is observed between the actual measurements ${{\mathrm {z}}}$ and the measurement ${{\mathrm {H}}}({\widehat{{{\mathrm {y}}}}})$ that would have been produced by ${\widehat{{{\mathrm {y}}}}}$. This discrepancy

$$\begin{aligned} {{\mathrm {J}}}({\widehat{{{\mathrm {y}}}}},t) = {{\mathrm {z}}}- {{\mathrm {H}}}({\widehat{{{\mathrm {y}}}}}) \end{aligned}$$

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 measurements—as 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 ${{\mathrm {F}}}({{\mathrm {y}}},t) = {{\mathrm {A}}}(t) {{\mathrm {y}}}+ {{\mathrm {R}}}$ and ${{\mathrm {H}}}({{\mathrm {y}}},t) = {{\mathrm {H}}}(t) {{\mathrm {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 ${{\mathrm {G}}}(t) = {{\mathrm {P}}}(t){{\mathrm {H}}}(t)^*$ where ${{\mathrm {P}}}$ is an operator following the Riccati evolution equation

$$\begin{aligned} \dot{{{\mathrm {P}}}} = {{\mathrm {A}}}{{\mathrm {P}}} + {{\mathrm {P}}}{{{\mathrm {A}}}}^{*} - {{\mathrm {P}}}{{\mathrm {H}}}^*{{\mathrm {H}}}{{\mathrm {P}}}, \quad {{\mathrm {P}}}(0) = {{\mathrm {P}}}_0, \end{aligned}$$

and ${{\mathrm {H}}}^*$ is the adjoint of ${{\mathrm {H}}}$. Although ${{\mathrm {P}}}$ is computable for any dynamics operator ${{\mathrm {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 ${\widetilde{{{\mathrm {y}}}}} = {{\mathrm {y}}}^{{\text {ref}}} - {\widehat{{{\mathrm {y}}}}}$, we get in a fully linear setting the following dynamics

$$\begin{aligned} {\left\{ \begin{array}{ll} \dot{{\widetilde{{{\mathrm {y}}}}}} = ({{\mathrm {A}}}- {{\mathrm {G}}} {{\mathrm {H}}}) {\widetilde{{{\mathrm {y}}}}} - {{\mathrm {G}}}\chi \\ {\widetilde{{{\mathrm {y}}}}}(0) = \zeta _{{\mathrm {y}}}. \end{array}\right. } \end{aligned}$$

Hence, ${{\mathrm {G}}}$ should be designed to stabilize the estimation error dynamics operator ${{\mathrm {A}}}- {{\mathrm {G}}} {{\mathrm {H}}}$, so that the homogeneous system tends to 0, namely,

$$\begin{aligned} {\widetilde{{{\mathrm {y}}}}}_{\chi =0} \xrightarrow {t \rightarrow +\infty } 0. \end{aligned}$$

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

$$\begin{aligned} {{\mathrm {G}}} = \gamma {{\mathrm {H}}}^*, \end{aligned}$$

(1)

with $\gamma $ 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 ${{\mathrm {H}}}^*$, by the adjoint of the tangent operator $D{{\mathrm {H}}}({\widehat{{{\mathrm {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 $\theta $ the uncertain parameter to be identified. Note that $\theta $ may be a vector of components or even a distributed field. The main idea is to introduce an augmented state vector and dynamics operator

$$\begin{aligned} {}^{\flat }{{\mathrm {y}}}= \begin{pmatrix} {{\mathrm {y}}}\\ \theta \end{pmatrix}, \quad {}^{\flat }{{\mathrm {F}}}({}^{\flat }{{\mathrm {y}}},t) = \begin{pmatrix} {{\mathrm {F}}}({{\mathrm {y}}},\theta ,t) \\ 0 \end{pmatrix}, \quad \end{aligned}$$

such that we still have ${}^{\flat }{\dot{{{\mathrm {y}}}}} = {}^{\flat }{{\mathrm {F}}}({}^{\flat }{{\mathrm {y}}},t)$. Then, a Kalman observer can be directly defined on this augmented model leading to a covariance operator and gain

$$\begin{aligned} {{\mathrm {P}}}= \begin{pmatrix} {{\mathrm {P}}}_{{{\mathrm {y}}}{{\mathrm {y}}}} &{} {{\mathrm {P}}}_{{{\mathrm {y}}}\theta } \\ {{\mathrm {P}}}_{\theta {{\mathrm {y}}}} &{} {{\mathrm {P}}}_{\theta \theta } \end{pmatrix}, \quad {{\mathrm {G}}}= \begin{pmatrix} {{\mathrm {G}}}_{{\mathrm {y}}}\\ {{\mathrm {G}}}_\theta \end{pmatrix}. \end{aligned}$$

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 $\gamma {{\mathrm {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

(2)

where ${{\mathrm {U}}}^{-1}$ is a reduced covariance operator on the parameter space and ${{\mathrm {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

$$\begin{aligned} {{\mathrm {G}}}_{{\mathrm {y}}}= (\gamma \mathbb {1} + {{\mathrm {L}}} {{\mathrm {U}}}^{-1} {{\mathrm {L}}}^*) D{{\mathrm {H}}}^*. \end{aligned}$$

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 $\Omega (t)$ at any time t. This domain is the image of a reference configuration $\Omega _0$ through the solid deformation mapping

$$\begin{aligned} \varvec{\varphi }: \left| \begin{aligned} \Omega _0\times [0,T]&\longrightarrow \Omega (t)\\ (\varvec{\xi },t)&\longmapsto \varvec{x}= \varvec{\varphi }(\varvec{\xi },t) = \varvec{\xi }+ \varvec{u}(\varvec{\xi },t), \end{aligned} \right. \end{aligned}$$

where $\varvec{u}$ denotes the solid displacement, so that the solid velocity is given by $\varvec{v}= {\dot{\varvec{u}}}$. The deformation gradient tensor $\varvec{F}$ is given by

$$\begin{aligned} \varvec{F}(\varvec{\xi },t) = \varvec{\nabla }_{\varvec{\xi }} \varvec{\varphi }= \varvec{\mathbb {1}} + \varvec{\nabla }_{\varvec{\xi }}\varvec{u}. \end{aligned}$$

Furthermore, we introduce the right Cauchy-Green deformation tensor $\varvec{C}= \varvec{F}^{\intercal }\cdot \varvec{F}$. We finally recall that the Green-Lagrange strain tensor denoted by $\varvec{e}$ is defined by

$$\begin{aligned} \varvec{e}= \frac{1}{2} (\varvec{C}-\varvec{\mathbb {1}}) = \frac{1}{2} \Bigl (\varvec{\nabla }_{\varvec{\xi }}\varvec{u}+ (\varvec{\nabla }_{\varvec{\xi }}\varvec{u})^\intercal + (\varvec{\nabla }_{\varvec{\xi }}\varvec{u})^\intercal \cdot \varvec{\nabla }_{\varvec{\xi }}\varvec{u}\Bigr ). \end{aligned}$$

Regarding the mechanical quantities notation, we denote by $\rho $ the tissue mass per unit volume, and by $\varvec{\sigma }$ 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 $\varvec{T}= J \varvec{\sigma }\cdot \varvec{F}^{-\intercal }$ and $\varvec{\Sigma }= \varvec{F}^{-1} \cdot \varvec{T}= J \varvec{F}^{-1} \cdot \varvec{\sigma }\cdot \varvec{F}^{-\intercal }$, where $J = \det {\varvec{F}}$.

The constitutive law can be considered as a nonlinear rheological combination of a passive part and an active part $\varvec{T}= \varvec{\Sigma }_p + \varvec{\Sigma }_a$. The passive part $\varvec{\Sigma }_p$ is described by a hyperelastic law of potential ${\mathcal {W}}$ and a viscous component chosen proportional to the strain rate ${\dot{\varvec{e}}}$

$$\begin{aligned} \varvec{\Sigma }_p(\varvec{e},{\dot{\varvec{e}}}) = \frac{\partial {\mathcal {W}}}{\partial \varvec{e}}(\varvec{e}) + \eta _s {\dot{\varvec{e}}}. \end{aligned}$$

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 in-vivo 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.

$$\begin{aligned} {\mathcal {W}}= C_0\exp \big (C_1(J_1 - 3)^2\big ) + C_2\exp \big (C_3(J_4 - 1)^2\big ) + \kappa (J - 1) - \kappa \ln (J), \end{aligned}$$

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 $\varvec{\tau }$, namely $ J_1 = J^{-\frac{2}{3}}{{\text {tr}}}(\varvec{C}),~J_4 = J^{-\frac{2}{3}}\varvec{\tau }\cdot \varvec{C}\cdot \varvec{\tau }, $ and $\kappa $ is the bulk coefficient.

For the active part $\varvec{\Sigma }_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 $\tau _c$, along the fiber direction $\varvec{\tau }$ in a chemically-controlled constitutive law describing myofibre mechanics [27, 28]. Therefore, we have $\varvec{\Sigma }_a = \sigma _{{\scriptscriptstyle 1{{\text {D}}}}}(e_c,k_c,\tau _c) \varvec{\tau }\otimes \varvec{\tau }$. We finally end up with the following second Piola-Kirchhoff stress tensor

$$\begin{aligned} \varvec{\Sigma }(\varvec{e},e_c,k_c,\tau _c) = \frac{\partial {\mathcal {W}}}{\partial \varvec{e}}(\varvec{e}) + \eta _s {\dot{\varvec{e}}} + \sigma _{{\scriptscriptstyle 1{{\text {D}}}}}(e_c,k_c,\tau _c)\varvec{\tau }\otimes \varvec{\tau }. \end{aligned}$$

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 $ \varvec{T}\cdot \varvec{n}= k_s \varvec{u}+ c_s \varvec{v}{{\text { on }}} \Gamma _n$, where $\varvec{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 $ \varvec{\sigma }\cdot \varvec{n}_t = - p_{v,i} \, \varvec{n}_t {{\text { on }}} \Gamma _{n,i}(t), i=\{1,2\}, $ where, here, $\varvec{n}_t$ denotes the normal of the deformed configuration boundary. Finally, the complete mechanical model reads

$$\begin{aligned} \left\{ \begin{array}{ll} {\dot{\varvec{u}}} = \varvec{v}, &{} \quad \text{ in } \Omega _0\\ \rho {\dot{\varvec{v}}} - \mathop {\mathbf {div}}(\varvec{T}) = 0,&{} \quad \text{ in } \Omega _0\\ \varvec{T}\cdot \varvec{n}= k_s \varvec{u}+ c_s \varvec{v}, &{} \quad {{\text {on }}} \Gamma _n \\ \varvec{T}\cdot \varvec{n}= - J p_{v,i} \, \varvec{F}^{-\intercal } \cdot \varvec{n}, &{} \quad {{\text {on }}} \Gamma _{c,i} \\ \varvec{T}\cdot \varvec{n}= 0, &{} \quad {{\text {on }}} \partial \Omega _0\backslash ((\cup _i \Gamma _{c,i}) \cup \Gamma _n). \end{array} \right. \end{aligned}$$

(3)

Together with the internal variable dynamics given in [27], it constitutes a general definition of the dynamical operator denoted by ${{\mathrm {F}}}$ in our above summarized description for a state ${{\mathrm {y}}}$ corresponding to $(\varvec{u},\varvec{v},e_c,k_c,\tau _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 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 $\omega _0$ of the geometry. We introduce the observation operator ${{\mathrm {H}}}= ({{\mathrm {H}}}_{|{\mathcal {Y}}^{u}} \quad 0) \in {\mathcal {L}}({\mathcal {Y}},{\mathcal {Z}})$ such that

$$\begin{aligned} {{\mathrm {H}}}_{|{\mathcal {Y}}^{u}} : \left| \begin{aligned} {\mathcal {Y}}^{u}&\rightarrow {\mathcal {Z}}\\ \varvec{u}&\mapsto \mathbb {1}_{\omega _0} \varvec{u}, \end{aligned} \right. \end{aligned}$$

where ${\mathcal {Y}}^{u}= {\mathcal {H}}^1(\Omega _0)^3$ is the displacement space. Note that ${{\mathrm {H}}}$ apply to the state space ${\mathcal {Y}}$ gathering the displacement space and the velocity space (and eventually additional variables in more general mechanical configuration). By contrast, ${{\mathrm {H}}}_{|{\mathcal {Y}}^{u}}$ is the same operator of input space restricted to the displacement space ${\mathcal {Y}}^{u}$ . We then need to specify the observation space ${\mathcal {Z}}$. One possible choice is to consider ${\mathcal {L}}^2(\omega _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 ${\mathcal {H}}^1(\Omega _0)^3$, the space of square-integrable vector fields with square-integrable first derivatives. We should rather consider ${\mathcal {Z}}= {\mathcal {H}}^1(\omega _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

$$\begin{aligned} {\left\{ \begin{array}{ll} - \mathop {\mathbf {div}}(\varvec{\sigma }^{{{\text {lin}}}}(\varvec{\psi })) = 0, &{}\quad {{\text {in}}} ~\Omega _0\backslash \omega _0\\ \varvec{\psi } = \varvec{\varphi }, &{} \quad {{\text {in}}} ~\omega _0\\ \varvec{\sigma }^{{{\text {lin}}}}(\varvec{\psi }) \cdot \varvec{n}= k_s \varvec{\psi }, &{}\quad {{\text {on }}} \Gamma _n \\ \varvec{\sigma }^{{{\text {lin}}}}(\varvec{\psi }) \cdot \varvec{n}= 0, &{}\quad {{\text {on }}} \partial \Omega _0\backslash \Gamma _n, \end{array}\right. } \end{aligned}$$

(4)

where $\varvec{\sigma }^{{{\text {lin}}}}$ denotes the stress tensor given by linearized isotropic elasticity. In particular, we denote the corresponding linear constitutive law by

$$\begin{aligned} \varvec{\sigma }^{{{\text {lin}}}}(\varvec{\psi }) = \varvec{A}: \varvec{\varepsilon }(\varvec{\psi }), \end{aligned}$$

where $\varvec{\varepsilon }$ denotes the usual linearized strain tensor. The boundary conditions on $\partial \Omega _0$ are adequately obtained from (3). In the following, we will denote the extension operator by

$$\begin{aligned} \varvec{\psi } = \mathbf {Ext}_{\omega _0}(\varvec{\varphi }). \end{aligned}$$

Note that an equivalent variational characterization of the extension is given by

$$\begin{aligned} \forall \varvec{v}^{\sharp }\in {\mathcal {Y}}^{u}\quad {{\text {s.t. }}} \varvec{v}^{\sharp }_{|\omega _0} = 0, \quad (\mathbf {Ext}_{\omega _0}(\varvec{\varphi }),\varvec{v}^{\sharp })_{{\mathcal {E}}_0}=0, \end{aligned}$$

(5)

where the energy dot-product is here defined by

$$\begin{aligned} (\varvec{u}_1,\varvec{u}_2)_{{\mathcal {E}}_0}&= \int _{\Omega _0} \varvec{\sigma }^{{{\text {lin}}}}(\varvec{u}_1) : \varvec{\varepsilon }(\varvec{u}_2) \, {{\mathrm {d}}}\Omega + \int _{\Gamma _n} k_s \, \varvec{u}_1\cdot \varvec{u}_2 \, {{\mathrm {d}}}\Gamma . \end{aligned}$$

We can prove—see Appendix A for a similar dot product $(\cdot ,\cdot )_{{\mathcal {E}}_0}$—that

$$\begin{aligned} (\mathbf {Ext}_{\omega _0}(\varvec{\varphi }),\mathbf {Ext}_{\omega _0}(\varvec{\varphi }))_{{\mathcal {E}}_0}^{\frac{1}{2}} \end{aligned}$$

is a norm in ${\mathcal {Z}}= {\mathcal {H}}^1(\omega _0)^3$. 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 ${{\mathrm {H}}}^*$ is given by

$$\begin{aligned} {{\mathrm {H}}}^* = \begin{pmatrix} {{\mathrm {H}}}_{|{\mathcal {Y}}^{u}}^* \\ 0 \end{pmatrix}, {{\text { with }}} {{\mathrm {H}}}_{|{\mathcal {Y}}^{u}}^* : \left| \begin{aligned} {\mathcal {Z}}&\rightarrow {\mathcal {Y}}^{u}\\ \varvec{{{\mathrm {z}}}}&\mapsto \mathbf {Ext}_{\omega _0}(\varvec{{{\mathrm {z}}}}). \end{aligned} \right. \end{aligned}$$

Note that in the rest of the document, by a slight abuse of notation, we will denote by ${{\mathrm {H}}}$ the operator applying either on the state space or on the displacement field extracted from the complete state ${{\mathrm {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 = \gamma {{\mathrm {H}}}^*$ is given in strong formulation by:

(6)

where $\Gamma _0 = \partial \Omega _0\backslash ((\cup _i \Gamma _{c,i}) \cup \Gamma _n)$. 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

$$\begin{aligned} {\widetilde{{{\mathrm {y}}}}} = \big ( {\widetilde{\varvec{u}}}, ~{\widetilde{\varvec{v}}} \big )^\intercal = \big (\varvec{u}^{{\text {ref}}} - {\widehat{\varvec{u}}}, ~\varvec{v}^{{\text {ref}}} - {\widehat{\varvec{v}}} \big )^\intercal , \end{aligned}$$

in the simplified case of linearized visco-elasticity without activation internal variables, namely

$$\begin{aligned} \varvec{\sigma }= \varvec{\sigma }^{{{\text {lin}}}}(\varvec{u}) + \eta _s \varvec{\varepsilon }(\varvec{v}). \end{aligned}$$

The estimation error satisfies the following weak formulation, for any $\varvec{v}^{\sharp }\in {\mathcal {Y}}^{u},$

$$\begin{aligned} \int _{\Omega _0} \rho \, \dot{{\widetilde{\varvec{v}}}}\cdot \varvec{v}^{\sharp }\, {{\mathrm {d}}}\Omega + ({\widetilde{\varvec{u}}},\varvec{v}^{\sharp })_{{\mathcal {E}}_0} + \int _{\Omega _0} \eta _s \, \varvec{\varepsilon }({\widetilde{\varvec{v}}}) : \varvec{\varepsilon }(\varvec{v}^{\sharp }) \, {{\mathrm {d}}}\Omega + \int _{\Gamma _n} c_s \, {\widetilde{\varvec{v}}}\cdot \varvec{v}^{\sharp }\, {{\mathrm {d}}}\Gamma = 0, \end{aligned}$$

with the additional observer-based relation

$$\begin{aligned} \dot{{\widetilde{\varvec{u}}}} = {\widetilde{\varvec{v}}} - \gamma \,\mathbf {Ext}_{\omega _0}(\mathbb {1}_{\omega _0} {\widetilde{\varvec{u}}}), \end{aligned}$$

assuming zero measurement error to fix the ideas. Weighing the latter relation by ${\widetilde{\varvec{u}}}$ and using the energy dot-product yields

$$\begin{aligned} ({\widetilde{\varvec{v}}},{\widetilde{\varvec{u}}})_{{\mathcal {E}}_0} = \frac{1}{2}\frac{{{\mathrm {d}}}}{{{\mathrm {d}}} t} \left\| {\widetilde{\varvec{u}}}\right\| _{{\mathcal {E}}_0}^2 + \gamma (\mathbf {Ext}_{\omega _0}(\mathbb {1}_{\omega _0}{\widetilde{\varvec{u}}}),{\widetilde{\varvec{u}}})_{{\mathcal {E}}_0} = \frac{1}{2}\frac{{{\mathrm {d}}}}{{{\mathrm {d}}} t} \left\| {\widetilde{\varvec{u}}}\right\| _{{\mathcal {E}}_0}^2 + \gamma \left\| \mathbf {Ext}_{\omega _0}(\mathbb {1}_{\omega _0}{\widetilde{\varvec{u}}})\right\| _{{\mathcal {E}}_0}^2, \end{aligned}$$

where we have used the orthogonality property (5). We can now substitute this expression in the above variational formulation applied with the test function $\varvec{v}^{\sharp }={\widetilde{\varvec{v}}}$, which gives

$$\begin{aligned}&\frac{{{\mathrm {d}}}}{{{\mathrm {d}}} t} \frac{1}{2}\bigl ( \left\| {\widetilde{\varvec{v}}}\right\| _{{\mathcal {K}}}^2 + \left\| {\widetilde{\varvec{u}}}\right\| _{{\mathcal {E}}_0}^2\bigr ) \nonumber \\&\quad =- \int _{\Omega _0} \eta _s \, \varvec{\varepsilon }({\widetilde{\varvec{v}}}) : \varvec{\varepsilon }({\widetilde{\varvec{v}}}) \, {{\mathrm {d}}}\Omega - \int _{\Gamma _n} c_s \, {\widetilde{\varvec{v}}}\,^2 \, {{\mathrm {d}}}\Gamma - \gamma \left\| \mathbf {Ext}_{\omega _0}(\mathbb {1}_{\omega _0}{\widetilde{\varvec{u}}})\right\| _{{\mathcal {E}}_0}^2, \end{aligned}$$

(7)

where $\left\| {\widetilde{\varvec{v}}}\right\| _{{\mathcal {K}}}^2 = (\rho {\widetilde{\varvec{v}}}, {\widetilde{\varvec{v}}})_{{\mathcal {L}}^2(\Omega _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 ${\mathcal {S}}_0$ a surface embedded in the reference domain $\Omega _0$ and $\varvec{e}$ a vector field defined on ${\mathcal {S}}_0$, with $(\varvec{e}_{\bot }^1,\varvec{e}_{\bot }^2)$ defined so that $(\varvec{e}_{\bot }^1,\varvec{e}_{\bot }^2,\varvec{e})$ gives an orthonormal basis at any point in ${\mathcal {S}}_0$. We define the extension $\varvec{\psi } = \mathbf {Ext}_{{\mathcal {S}}_0}(\varvec{e}~;~\varphi )$ from ${\mathcal {S}}_0$ of the scalar field $\varphi $ in the direction $\varvec{e}$ as

$$\begin{aligned} {\left\{ \begin{array}{ll} - \mathop {\mathbf {div}}(\varvec{\sigma }^{{{\text {lin}}}}(\varvec{\psi })) = 0, &{} \quad {{\text {in}}} ~\Omega _0\\ \varvec{\psi }\cdot \varvec{e}= \varphi , &{} \quad {{\text {on}}} ~{\mathcal {S}}_0\\ (\varvec{\sigma }^{{{\text {lin}}}}(\varvec{\psi })\cdot \varvec{n})\cdot \varvec{e}_{\bot }^1 = (\varvec{\sigma }^{{{\text {lin}}}}\cdot \varvec{n})\cdot \varvec{e}_{\bot }^2= 0, &{}\quad {{\text {on}}}~{\mathcal {S}}_0\\ \varvec{\sigma }^{{{\text {lin}}}}(\varvec{\psi }) \cdot \varvec{n}= k_s \varvec{\psi }, &{} \quad {{\text {on }}} \Gamma _n \\ \varvec{\sigma }^{{{\text {lin}}}}(\varvec{\psi }) \cdot \varvec{n}= 0, &{} \quad {{\text {on }}} \partial \Omega _0\backslash \Gamma _n. \end{array}\right. } \end{aligned}$$

(8)

An equivalent variational characterization of this extension operator is

$$\begin{aligned} \forall \varvec{v}^{\sharp }\in {\mathcal {Y}}^{u}\quad {{\text {s.t. }}} \varvec{v}^{\sharp }\cdot \varvec{e}=0 {{\text { on }}} {\mathcal {S}}_0, \quad (\mathbf {Ext}_{{\mathcal {S}}_0}(\varvec{e}~;~\varphi ),\varvec{v}^{\sharp })_{{\mathcal {E}}_0}=0. \end{aligned}$$

(9)

We can define a norm on the surface-based data using this extension, namely, $\left\| \varvec{\psi }\right\| _{{\mathcal {E}}_0}$. Typically, considering the following linear observation operator

$$\begin{aligned} {{\mathrm {H}}}: \left| \begin{aligned} {\mathcal {Y}}^{u}&\rightarrow {\mathcal {Z}}\\ \varvec{u}&\mapsto \varvec{u}|_{{\mathcal {S}}_0}\cdot \varvec{e}, \end{aligned} \right. \end{aligned}$$

we can use this norm in the observation space. Accordingly, observer terms in the form ${{\mathrm {H}}}^*({{\mathrm {z}}}-{{\mathrm {H}}}\, \widehat{\varvec{u}})$ will give in a variational setting

$$\begin{aligned} (\mathbf {Ext}_{{\mathcal {S}}_0}(\varvec{e}~;~{{\mathrm {z}}}-{{\mathrm {H}}}\, \widehat{\varvec{u}}),\mathbf {Ext}_{{\mathcal {S}}_0}(\varvec{e}~;~{{\mathrm {H}}}\, \varvec{v}^{\sharp }))_{{\mathcal {E}}_0} = (\mathbf {Ext}_{{\mathcal {S}}_0}(\varvec{e}~;~{{\mathrm {z}}}-{{\mathrm {H}}}\, \widehat{\varvec{u}}),\varvec{v}^{\sharp })_{{\mathcal {E}}_0}, \end{aligned}$$

where the second expression is obtained by using the characterization (9) when observing that on ${\mathcal {S}}_0$

$$\begin{aligned} \mathbf {Ext}_{{\mathcal {S}}_0}(\varvec{e}~;~{{\mathrm {H}}}\, \varvec{v}^{\sharp }) \cdot \varvec{e} = \varvec{v}^{\sharp }\cdot \varvec{e}. \end{aligned}$$

Therefore, we have in this case

$$\begin{aligned} {{\mathrm {H}}}^*({{\mathrm {z}}}-{{\mathrm {H}}}\, \widehat{\varvec{u}}) = \mathbf {Ext}_{{\mathcal {S}}_0}(\varvec{e}~;~{{\mathrm {z}}}-{{\mathrm {H}}}\, \widehat{\varvec{u}}). \end{aligned}$$

(10)

Now, when dealing more generally with a discrepancy operator ${{\mathrm {J}}}(\widehat{\varvec{u}},{{\mathrm {z}}})$ pertaining to the displacements on a surface ${\mathcal {S}}_0$, we will generalize this strategy by using the observer correction given by

$$\begin{aligned} -\mathbf {Ext}_{{\mathcal {S}}_0}(D_u{{\mathrm {J}}}(\widehat{\varvec{u}},{{\mathrm {z}}})~;~{{\mathrm {J}}}(\widehat{\varvec{u}},{{\mathrm {z}}})), \end{aligned}$$

obtained by directly substituting in (10) ${{\mathrm {J}}}(\widehat{\varvec{u}},{{\mathrm {z}}})$ for $({{\mathrm {z}}}-{{\mathrm {H}}}\, \widehat{\varvec{u}})$, and $-D_u{{\mathrm {J}}}(\widehat{\varvec{u}},{{\mathrm {z}}})$—associated with a vector field on ${\mathcal {S}}_0$—for $\varvec{e}$ since in the linear case we have ${{\mathrm {H}}}\, \widehat{\varvec{u}}=\widehat{\varvec{u}}\cdot \varvec{e}$ on ${\mathcal {S}}_0$. Of course, this also easily generalizes to a measurement made of a collection of such surface-based data associated with $N_{\mathcal {S}}$ surfaces $({\mathcal {S}}_0^i)_{i=1}^{N_{\mathcal {S}}}$ and associated vector fields $\varvec{e}^i$, for which the correction will be given (in the linear case) by

$$\begin{aligned} {{\mathrm {H}}}^*({{\mathrm {z}}}-{{\mathrm {H}}}\,\widehat{\varvec{u}}) = \sum _{i=1}^{N_{\mathcal {S}}}\mathbf {Ext}_{{\mathcal {S}}_0^i} (\varvec{e}^i~;~{{\mathrm {z}}}-{{\mathrm {H}}}\,\widehat{\varvec{u}}). \end{aligned}$$

(11)

Tag planes

We consider data consisting of a set of $N_{{\mathcal {P}}}$ tag planes ${\mathcal {T}}= \mathop {\bigcup }\limits _i^{N_{{\mathcal {P}}}}{\mathcal {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 ${\mathcal {T}}_0 = \mathop {\bigcup }\limits _i^{N_{{\mathcal {P}}}}{\mathcal {P}}_0^i$ the set of tag planes in the reference configuration, mapped by the estimated trajectory $\widehat{\varvec{u}}$ to ${\widehat{{\mathcal {T}}}}= \mathop {\bigcup }\limits _i^{N_{{\mathcal {P}}}}{\widehat{{\mathcal {P}}}}^i$. For any point in a model tag plane $\widehat{\varvec{x}}= \varvec{\xi }+\widehat{\varvec{u}}(\varvec{\xi })\in {\widehat{{\mathcal {P}}}}^i$ for some $\varvec{\xi }\in {\mathcal {P}}_0^i$, we can compute the signed distance to the corresponding target tag plane ${\mathcal {P}}^i$ by

$$\begin{aligned} {{\mathrm {dist}}}(\widehat{\varvec{x}},{\mathcal {P}}^i) = (\widehat{\varvec{x}}- \Pi _{{\mathcal {P}}^i}\widehat{\varvec{x}})\cdot \varvec{n}_{{\mathcal {P}}^i}, \end{aligned}$$

(12)

where $\widehat{\varvec{x}}$ is a point on the ith model tag plane ${\widehat{{\mathcal {P}}}}^i$, $\Pi _{{\mathcal {P}}^i}\widehat{\varvec{x}}$ is the projection of this point on the corresponding observed tag plane, and $\varvec{n}_{{\mathcal {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 ${\mathcal {T}}_0$ . When differentiating with respect to the displacement field we have

$$\begin{aligned} \forall \varvec{v}^{\sharp }\in {\mathcal {Y}}^{u},\quad D_u{{\mathrm {dist}}}(\widehat{\varvec{x}},{\mathcal {P}}^i)\cdot \varvec{v}^{\sharp }= \varvec{n}_{{\mathcal {P}}^i} \cdot \varvec{v}^{\sharp }. \end{aligned}$$

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 into

(13)

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_{{\mathcal {P}}}$ such lines $\bigl ({\mathcal {L}}^{ij}\bigr )_{i=1}^{N_{{\mathcal {P}}}}$ in each 2D image ${\mathcal {I}}^j$, with $1\le j \le N_{{\mathcal {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 $\Pi _{{\mathcal {I}}^j}$ the Euclidean orthogonal projection onto the image ${\mathcal {I}}^j$, we compute the signed distance of the projected point to the corresponding tag line within each image—see Fig. 4—i.e.

$$\begin{aligned} {{\mathrm {dist}}}(\Pi _{{\mathcal {I}}^j}\widehat{\varvec{x}},{\mathcal {L}}^{ij}) = \bigl ( \Pi _{{\mathcal {I}}^j}\widehat{\varvec{x}}- \Pi _{{\mathcal {L}}^{ij}}\Pi _{{\mathcal {I}}^j}\widehat{\varvec{x}}\bigr ) \cdot \varvec{n}_{{\mathcal {L}}^{ij}}. \end{aligned}$$

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,

$$\begin{aligned} {{\mathrm {J}}}_i(\widehat{\varvec{u}},{{\mathrm {z}}}) = {\mathcal {J}}_{(j)}\bigl ({{\mathrm {dist}}}(\Pi _{{\mathcal {I}}^j}\widehat{\varvec{x}},{\mathcal {L}}^{ij})\bigr ), \end{aligned}$$

(14)

for each plane ${\mathcal {P}}_0^i$, where ${\mathcal {J}}_{(j)}$ denotes the interpolation operator. When differentiating this expression, we have

$$\begin{aligned} \forall \varvec{v}^{\sharp }\in {\mathcal {Y}}^{u},\quad D_u{{\mathrm {dist}}}(\Pi _{{\mathcal {I}}^j}\widehat{\varvec{x}},{\mathcal {L}}^{ij})\cdot \varvec{v}^{\sharp }= \varvec{n}_{{\mathcal {L}}^{ij}} \cdot \varvec{v}^{\sharp }, \end{aligned}$$

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 ${\mathcal {J}}_{(j)}'$ the derivative with respect to this coordinate, a straightforward computation finally yields

$$\begin{aligned} \forall \varvec{v}^{\sharp }\in {\mathcal {Y}}^{u},\quad D_u{{\mathrm {J}}}_i(\widehat{\varvec{u}},{{\mathrm {z}}}) \cdot \varvec{v}^{\sharp }= \Bigl ( {\mathcal {J}}_{(j)}\bigl (\varvec{n}_{{\mathcal {L}}^{ij}}\bigr ) + {\mathcal {J}}_{(j)}' \bigl ({{\mathrm {dist}}}(\Pi _{{\mathcal {I}}^j}\widehat{\varvec{x}},{\mathcal {L}}^{ij})\bigr ) \varvec{n}_{{\mathcal {I}}} \Bigr ) \cdot \varvec{v}^{\sharp }, \end{aligned}$$

where $\varvec{n}_{{\mathcal {I}}}$ denotes the direction vector of the orthogonal axis. Note that when considering e.g. linear interpolation, the derivative ${\mathcal {J}}_{(j)}'$ is directly given by the finite difference expression computed between the two adjacent planes. This finally gives for the observer correction equation

(15)

with $\varvec{e}^i={\mathcal {J}}_{(j)}\bigl (\varvec{n}_{{\mathcal {L}}^{ij}}\bigr ) + {\mathcal {J}}_{(j)}' \bigl ({{\mathrm {dist}}}(\Pi _{{\mathcal {I}}^j}\widehat{\varvec{x}},{\mathcal {L}}^{ij})\bigr ) \,\varvec{n}_{{\mathcal {I}}}$.

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-the-plane 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 $D{{\mathrm {H}}}$, 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-component field over a plane—instead of a scalar field for the above distances. Hence, we resort to a slightly different extension operator, namely, $\varvec{\psi } = \mathbf {Ext}_{{\mathcal {S}}_0}'(\varvec{e}~;~\varvec{\varphi })$ defined by

$$\begin{aligned} {\left\{ \begin{array}{ll} - \mathop {\mathbf {div}}(\varvec{\sigma }^{{{\text {lin}}}}(\varvec{\psi })) = 0, &{} \quad {{\text {in}}} ~\Omega _0\\ \Pi _{\varvec{e}}\varvec{\psi } = \varvec{\varphi }, &{} \quad {{\text {on}}} ~{\mathcal {S}}_0\\ (\varvec{\sigma }^{{{\text {lin}}}}\cdot \varvec{n})\cdot \varvec{e} = 0, &{}\quad {{\text {on}}}~{\mathcal {S}}_0\\ \varvec{\sigma }^{{{\text {lin}}}} \cdot \varvec{n}= k_s \varvec{\psi }, &{}\quad {{\text {on }}} \Gamma _n \\ \varvec{\sigma }^{{{\text {lin}}}} \cdot \varvec{n}= 0, &{} \quad {{\text {on }}} \partial \Omega _0\backslash \Gamma _n, \end{array}\right. } \end{aligned}$$

(16)

where $\Pi _{\varvec{e}}$ denotes the projection onto the plane orthogonal to $\varvec{e}$, plane in which $\varvec{\varphi }$ is assumed to lie. Finally, the correction equation for the observer reads

(17)

where $\varvec{n}_{\mathcal {I}}$ denotes the normal to the image planes, and the innovation term ${{\mathrm {z}}}-{{\mathrm {H}}}(\widehat{\varvec{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 using a prediction–correction scheme

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 part—i.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

$$\begin{aligned} {{\mathrm {z}}}^n = {{\mathrm {H}}}({{\mathrm {y}}}(n\Delta {{\mathrm {t}}})), \end{aligned}$$

(18)

where ${{\mathrm {H}}}(\cdot )$ is a nonlinear and sufficiently smooth observation operator. Following [6] we propose to define the time-discrete observer as

where we denote by ${\widehat{{{\mathrm {y}}}}}^{n+1}_-$ and ${\widehat{{{\mathrm {y}}}}}^{n+1}_+$ the prediction and correction steps respectively. In the correction relation (19b), we use ${\widehat{{{\mathrm {y}}}}}^{e}_+$ an extrapolated trajectory, and $D{{\mathrm {H}}}^{e}_+$ the tangent operator of the observation operator evaluated at ${\widehat{{{\mathrm {y}}}}}^{e}_+$, i.e. $D{{\mathrm {H}}}^{e}_+ = D{{\mathrm {H}}}({\widehat{{{\mathrm {y}}}}}^{e}_+)$. While other choices are possible—see [6]—the most simple approach is to set ${\widehat{{{\mathrm {y}}}}}^{e}_+= {\widehat{{{\mathrm {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

$$\begin{aligned} {{\mathrm {J}}}\big ({{\mathrm {y}}}(n\Delta {{\mathrm {t}}}),{{\mathrm {z}}}^n\big ) = 0. \end{aligned}$$

In this case, the time-discrete observer can be directly inferred from system (19a)–(19b)

The corresponding linearized estimation error satisfies the estimate in Proposition B.4 of Appendix B, but the operator $D{{\mathrm {J}}}^{e}_+ = D{{\mathrm {J}}}({\widehat{{{\mathrm {y}}}}}^e_+,{{\mathrm {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\in {\mathbb {R}}^{N_{{\text {dof}}}}$ the vector of degrees of freedom—of dimension $N_{{{\text {dof}}}}$—associated with the function $\varvec{u}_h$ defined in the finite-element space ${\mathcal {Y}}^{u}_h \subset {\mathcal {Y}}^{u}$. We also define a state vector as the concatenation of the displacement degrees of freedom and the velocity degrees of freedom $Y= \bigl ( {\begin{matrix} U\\ V\end{matrix}} \bigr )$. This final space discretization procedure enables us to define a fully-discrete transition operator $F_{n+1|n}$, such that the prediction step (19a)—or (20a)—can be rewritten as

$$\begin{aligned} {\widehat{Y}}^{n+1}_- = F_{n+1|n}({\widehat{Y}}^{n}_+), \end{aligned}$$

(21)

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 $(\varvec{u}_{h,1}, \varvec{u}_{h,2}) \in {\mathcal {Y}}^{u}_h \times {\mathcal {Y}}^{u}_h$

$$\begin{aligned} U_1^\intercal MU_2 = (\rho \varvec{u}_{h,1},\varvec{u}_{h,2})_{{\mathcal {L}}^2(\Omega _0)^3} = \int _{\Omega _0} \rho \varvec{u}_{h,1} \cdot \varvec{u}_{h,2} \, {{\mathrm {d}}}\Omega , \end{aligned}$$

and

$$\begin{aligned} U_1^\intercal KU_2 = (\varvec{u}_{h,1},\varvec{u}_{h,2})_{{\mathcal {E}}_0} = \int _{\Omega _0} \varvec{\varepsilon }(\varvec{u}_{h,1}) : \varvec{A}: \varvec{\varepsilon }(\varvec{u}_{h,2}) \, {{\mathrm {d}}}\Omega + \int _{\Gamma _n} k_s \, \varvec{u}_{h,1} \cdot \varvec{u}_{h,2} \, {{\mathrm {d}}}\Gamma , \end{aligned}$$

the mass and stiffness matrix respectively, and we consider the associated norm

$$\begin{aligned} N= \begin{pmatrix} K&{} 0 \\ 0 &{} M\end{pmatrix}. \end{aligned}$$

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 $\bigl (\!{\begin{matrix} H&0 \end{matrix}} \!\bigr )$—or $\bigl (\!{\begin{matrix} DH&0 \end{matrix}} \!\bigr )$.

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 ${\mathcal {Z}}_{h^\prime }\subset {\mathcal {Z}}$ a suitable discretization of the observation space. The (linear) observation operator boils down to an interpolator between the mesh and the observation subdomain $\omega _0$, whereas S is defined for all $(\psi _{h^\prime ,1}, \psi _{h^\prime ,2}) \in {\mathcal {Z}}_{h^\prime }\times {\mathcal {Z}}_{h^\prime }$ by

$$\begin{aligned} \Psi _{1}^\intercal S\Psi _{2} = (\psi _{h^\prime ,1},\psi _{h^\prime ,1})_{{\mathcal {Z}}} = (\mathbf {Ext}_{\omega _0}(\psi _{h^\prime ,1}),\mathbf {Ext}_{\omega _0} (\psi _{h^\prime ,2}))_{{\mathcal {E}}_0}. \end{aligned}$$

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

$$\begin{aligned} N\Bigl ( \dfrac{{\widehat{Y}}^{n+1}_+ - {\widehat{Y}}^{n+1}_-}{\Delta {{\mathrm {t}}}} \Bigr ) = \gamma D{H^e_+}^\intercal S\big (Z^{n+1} - H({\widehat{Y}}^{e}_+) - D{H^e_+}({\widehat{Y}}^{n+1}_+ - {\widehat{Y}}^{e}_+) \big ), \end{aligned}$$

where $D{H^e_+}$ is the tangent of the observation operator around the extrapolated trajectory ${\widehat{Y}}^{e}_+$. To simplify the presentation, we set ${\widehat{Y}}^{e}_+~=~{\widehat{Y}}^{n+1}_-$, so that the correction step can be expressed in a more practical form as

$$\begin{aligned} {\widehat{Y}}^{n+1}_+ = {\widehat{Y}}^{n+1}_- + \gamma (N+ \gamma D{H^e_+}^\intercal S_{\Delta {{\mathrm {t}}}}D{H^e_+})^{-1} D{H^e_+}^\intercal S_{\Delta {{\mathrm {t}}}}(Z^{n+1} - H({\widehat{Y}}^{n+1}_-)), \end{aligned}$$

(22)

where $S_{\Delta {{\mathrm {t}}}} = \Delta {{\mathrm {t}}}S$. 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_{\Delta t}^{-1}$, the a priori state ${\widehat{Y}}^{n+1}_-$ and the a priori state covariance $\gamma N^{-1}$ [15]. Alternatively, it can be interpreted as solving the following minimization problem

$$\begin{aligned} \underset{{\widehat{Y}}}{{{\mathrm {min}}}}~\dfrac{1}{2}\Big \{~\Vert {\widehat{Y}}- {\widehat{Y}}^{n+1}_-\Vert ^2_N + \gamma \Vert Z^{n+1} - H({\widehat{Y}}^{n+1}_-) - D{H^e_+}({\widehat{Y}}- {\widehat{Y}}^{n+1}_-) \Vert ^2_{S_{\Delta {{\mathrm {t}}}}}\Big \}, \end{aligned}$$

where $\Vert \cdot \Vert _N$ and $\Vert \cdot \Vert _{S_{\Delta t}}$ are the norms induced by the matrices N and $S_{\Delta {{\mathrm {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 ${\widehat{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 $\gamma 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

$$\begin{aligned} N\le N+ \gamma D{H^e_+}^\intercal S_{\Delta {{\mathrm {t}}}}D{H^e_+} \le N\big (1 + O(\Delta {{\mathrm {t}}}) \big ), \end{aligned}$$

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 $D{H^e_+}^\intercal S_{\Delta {{\mathrm {t}}}}D{H^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+ \gamma D{H^e_+}^\intercal S_{\Delta {{\mathrm {t}}}}D{H^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 $D{H^e_+}^\intercal S_{\Delta {{\mathrm {t}}}}D{H^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 $\varvec{\psi } = {{\mathrm {H}}}^*\varvec{\varphi } = \mathbf {Ext}_{\omega _0}(\varvec{\varphi })$ for $\varvec{\varphi }\in {\mathcal {Z}}$ can be approximated, after space discretization, by solving

$$\begin{aligned} \min _{ \Psi \in {\mathbb {R}}^{N_{{\text {dof}}}} } \dfrac{1}{2}\Big \{~\varepsilon \Psi ^\intercal K\Psi + (\Phi - H\Psi )^\intercal M_{{{\text {obs}}}}(\Phi - H\Psi ) \Big \} , \end{aligned}$$

(23)

where $\varepsilon > 0$ is the (small) penalization parameter, $\Psi $ and $\Phi $ are vectors of degrees of freedom, $H$ is an interpolation operator between the mesh and the subdomain $\omega _0$, and $M_{{{\text {obs}}}}$ is the matrix associated with the ${\mathcal {L}}^2$-norm on the observation space

$$\begin{aligned} \Phi _{1}^\intercal M_{{{\text {obs}}}}\Phi _{2} = \int _{\omega _0} \varvec{\phi }_{h^\prime ,1} \cdot \varvec{\phi }_{h^\prime ,2} \, {{\mathrm {d}}}\Omega ,\quad \forall (\varvec{\phi }_{h^\prime ,1}, \varvec{\phi }_{h^\prime ,2})\in {\mathcal {Z}}_{h^\prime }\times {\mathcal {Z}}_{h^\prime }. \end{aligned}$$

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 $(\varvec{u}_{h,1}, \varvec{u}_{h,2}) \in {\mathcal {Y}}^{u}_h \times {\mathcal {Y}}^{u}_h$,

$$\begin{aligned} \begin{aligned} U_1^\intercal H^\intercal S HU_2&= (\mathbb {1}_{\omega _0}\varvec{u}_{h,1}, \mathbb {1}_{\omega _0}\varvec{u}_{h,2})_{{\mathcal {Z}}} \\&= (\mathbf {Ext}_{\omega _0}(\mathbb {1}_{\omega _0}\varvec{u}_{h,1}), \mathbf {Ext}_{\omega _0}(\mathbb {1}_{\omega _0}\varvec{u}_{h,2}))_{{\mathcal {E}}_0}\\&= (\varvec{u}_{h,1},\mathbf {Ext}_{\omega _0}(\mathbb {1}_{\omega _0}\varvec{u}_{h,2}))_{{\mathcal {E}}_0}, \end{aligned} \end{aligned}$$

where we have used the orthogonality property (5) to obtain the last relation. Then, replacing the extension operator by its penalized approximation (23) yields

$$\begin{aligned} U_1^\intercal H^\intercal S HU_2 = U_1^\intercal K(\varepsilon K+ H^\intercal M_{{{\text {obs}}}}H)^{-1} H^\intercal M_{{{\text {obs}}}}HU_2, \end{aligned}$$

from which we retrieve the relation given in [6],

$$\begin{aligned} H^\intercal S = K(\varepsilon K+ H^\intercal M_{{{\text {obs}}}}H)^{-1} H^\intercal M_{{{\text {obs}}}}. \end{aligned}$$

(24)

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 $\theta $ discretized in ${\mathbb {R}}^r$, we introduce $r+1$ so-called unitary simplex sampling points $I^{[i]}$ in the space ${\mathbb {R}}^{r}$ and the associated weights $\alpha _i$ with the following rules

$$\begin{aligned} \sum _{i=1}^{r+1} \alpha _i I_{[i]} = 0, \quad \sum _{i=1}^{r+1} \alpha _i I_{[i]} {I_{[i]}}^\intercal = \mathbb {1}, \end{aligned}$$

(25)

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_\alpha $ and precompute these unitary sigma-points. We also denote by $M_{\Delta t} = \Delta t M_{{{\text {obs}}}}$, and by $[I_{[*]}]$ the matrix concatenating the $(I_{[i]})_{1\le i\le r+1}$ vectors side by side, and similarly for other matrices aggregating some particle vectors.

We then perform at each time step

1.
Sampling ($1\le i \le r+1$):
$$\begin{aligned} {\left\{ \begin{array}{ll} Q_n = \sqrt{(U^n)^{-1}} \\ {\widehat{Y}}_{[i]+}^{n} = {\widehat{Y}}_{n}^+ + L^n_{{{\mathrm {y}}}} Q_n^\intercal I_{[i]} \\ {\widehat{\theta }}_{[i]+}^{n} = {\widehat{\theta }}_{n}^+ + L^n_{\theta } Q_n^\intercal I_{[i]} \end{array}\right. } \end{aligned}$$
(26a)
2.
State prediction ($1\le i \le r+1$):
$$\begin{aligned} {\left\{ \begin{array}{ll} {\widehat{Y}}_{[i]-}^{n+1} = F_{n+1|n}({\widehat{Y}}_{[i]+}^{n},{\widehat{\theta }}^{n}_{[i]+})\\ {\widehat{\theta }}_{[i]-}^{n+1} = {\widehat{\theta }}^{n}_{[i]+}\\ {\widehat{\theta }}_{-}^{n+1} = \sum \nolimits _{i=1}^{r+1} \alpha _i {\widehat{\theta }}^{n+1}_{[i]-} \end{array}\right. } \end{aligned}$$
(26b)
3.
State correction ($1\le i \le r+1$):
$$\begin{aligned} {\left\{ \begin{array}{ll} Z^{n+1}_{[i]-} = H^{n+1}({\widehat{Y}}_{[i]-}^{n+1}) \\ {\widehat{Y}}^{n+1}_{[i]-+} = {\widehat{Y}}^{n+1}_{[i]-} \\ \qquad \qquad + \gamma (N+ \gamma D{H^e_+}^\intercal S_{\Delta {{\mathrm {t}}}}D{H^e_+})^{-1} D{H^e_+}^\intercal S_{\Delta {{\mathrm {t}}}}(Z^{n+1} - Z^{n+1}_{[i]-}) \\ {\widehat{Y}}_{-+}^{n+1} = \sum \nolimits _{i=1}^{r+1} \alpha _i {\widehat{Y}}^{n+1}_{[i]-+} \end{array}\right. } \end{aligned}$$
(26c)
4.
Parametric correction:
$$\begin{aligned} {\left\{ \begin{array}{ll} L^{n+1}_{Y} = [{\widehat{Y}}_{[*]-+}^{n+1}] M_\alpha [I_{[*]}]^\intercal \\ L^{n+1}_{\theta } = [{\widehat{\theta }}_{[*]-+}^{n+1}] M_\alpha [I_{[*]}]^\intercal \\ Z^{n+1}_{-} = \sum \nolimits _{i=1}^{r+1} \alpha _i Z^{n+1}_{[i]-} \\ \Gamma ^{n+1} = [Z_{[*]-}^{n+1}] M_\alpha [I_{[*]}]^\intercal \\ U^{n+1} = \mathbb {1} + (\Gamma ^{n+1})^\intercal M_{\Delta t} \Gamma ^{n+1} \\ \Upsilon ^{n+1} = U^{n+1} (\Gamma ^{n+1})^\intercal M_{\Delta t} (Z^{n+1} - Z^{n+1}_{-}) \\ {\widehat{Y}}^{n+1}_{+} = {\widehat{Y}}^{n+1}_{-+} - L^{n+1}_{Y} \Upsilon ^{n+1} \\ {\widehat{\theta }}^{n+1}_{+} = {\widehat{\theta }}^{n+1}_{-} - L^{n+1}_{\theta } \Upsilon ^{n+1} \end{array}\right. } \end{aligned}$$
(26d)

Results

In order to illustrate and assess the above data assimilation method, we propose to perform parameter estimation in a synthetic data context. More precisely, in this applicative example we will extract from the direct simulation of an infarcted heart the tag planes, tag grids and apparent displacements. The infarct is represented in the model by increasing the stiffness and decreasing the contractility in the septum, see Fig. 6. To that purpose we define two parameters $\theta ^K$ and $\theta $ such that the constant values $(C_0,C_2)$, appearing in the hyperelastic potential, and the contractility of the tissue $\sigma _0$, appearing in the model proposed in [27], are transformed into

$$\begin{aligned} (C_0,C_2) \rightarrow 2^{\theta ^K} (C_0,C_2) \quad ,~\sigma _0 \rightarrow 2^{\theta }\sigma _0. \end{aligned}$$

This exponential dependence allows to keep these parameters positive, which is crucial from a physical standpoint, and to ensure stability of the numerical scheme during the estimation procedure.

Table 1 Set parameters for the biomechanical heart model (SI units)

Full size table

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 $(\theta ^K,\theta ) = (1,-1)$ in the septum, and $(\theta ^K,\theta ) = (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 ${\mathcal {J}}_{{\mathcal {P}}_{{{\mathrm {m}}}}^i}$ the interpolation operator of a single tag plane, its displacement $\varvec{u}_{{\mathcal {P}}_{{{\mathrm {m}}}}^i}$ is given in this context by

$$\begin{aligned} \varvec{u}_{{\mathcal {P}}_{{{\mathrm {m}}}}^i}(\varvec{\xi }_{{\mathrm {m}}},t) = \left\{ \begin{array}{ll} {\mathcal {J}}_{{\mathcal {P}}_{{{\mathrm {m}}}}^i}(\varvec{u}(\varvec{\xi },t) - \varvec{u}_{{{\mathrm {m}}}}), &{} \quad \text{ if } \varvec{\xi }_{{\mathrm {m}}}\in \Omega _{{\mathrm {m}}}\cap {\mathcal {P}}_{{{\mathrm {m}}}}^i \\ 0 , &{} \quad \text{ otherwise }, \end{array} \right. \end{aligned}$$

where $\varvec{u}_{{{\mathrm {m}}}}$ is the model displacement at marking time $t_{{\mathrm {m}}}$, $\varvec{\xi }_{{\mathrm {m}}}$ denotes the associated deformed coordinate, namely, $\varvec{\xi }_{{\mathrm {m}}}=\varvec{\varphi }(\varvec{\xi },t_{{\mathrm {m}}})$, and $\Omega _{{\mathrm {m}}}=\varvec{\varphi }(\Omega _0,t_{{\mathrm {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 regularization—we consider a set of elastic (3D) tag layers $\mathop {\bigcup }\limits _i^{N_{{\mathcal {P}}}}{\mathcal {V}}_{{{\mathrm {m}}}}^i$. In practice, each tag layer is built so that ${\mathcal {P}}_{{{\mathrm {m}}}}^i \subset {\mathcal {V}}_{{{\mathrm {m}}}}^i$. Hence, the displacement of a tag plane ${\mathcal {P}}_{{{\mathrm {m}}}}^i$ is derived from

$$\begin{aligned} \varvec{u}_{{\mathcal {P}}_{{{\mathrm {m}}}}^i}(\varvec{\xi }_{{\mathrm {m}}},t) =\varvec{u}_{{\mathcal {V}}_{{{\mathrm {m}}}}^i}(\varvec{\xi }_{{\mathrm {m}}},t) \vert _{{\mathcal {P}}_{{{\mathrm {m}}}}^i}, \end{aligned}$$

(27)

where $\varvec{u}_{{\mathcal {V}}_{{{\mathrm {m}}}}^i}$ is the displacement of the tag layer ${\mathcal {V}}_{{{\mathrm {m}}}}^i$, verifying

$$\begin{aligned} \left\{ \begin{array}{ll} - \mathop {\mathbf {div}}(\varvec{\sigma }^{{\mathcal {V}}}(\varvec{u}_{{\mathcal {V}}_{{{\mathrm {m}}}}^i})) = 0, &{}\quad \text{ in } {\mathcal {V}}_{{{\mathrm {m}}}}^i \backslash \big ( \Omega _{{\mathrm {m}}}\cap {\mathcal {P}}_{{{\mathrm {m}}}}^i \big )\\ \varvec{u}_{{\mathcal {V}}_{{{\mathrm {m}}}}^i} = {\mathcal {J}}_{{\mathcal {P}}_{{{\mathrm {m}}}}^i}(\varvec{u}- \varvec{u}_{{{\mathrm {m}}}}), &{} \quad \text{ in } \Omega _{{\mathrm {m}}}\cap {\mathcal {P}}_{{{\mathrm {m}}}}^i \\ \varvec{\sigma }^{{\mathcal {V}}}(\varvec{u}_{{\mathcal {V}}_{{{\mathrm {m}}}}^i})\cdot \varvec{n}= 0, &{} \quad \text{ on } \partial {\mathcal {V}}_{{{\mathrm {m}}}}^i \end{array} \right. \end{aligned}$$

(28)

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,

$$\begin{aligned} \varvec{u}_{{\mathcal {V}}_{{{\mathrm {m}}}}^i} = \mathop {{\mathbf {Ext}}}_{\Omega _{{\mathrm {m}}}\cap {\mathcal {P}}_{{{\mathrm {m}}}}^i}\big ({\mathcal {J}}_{{\mathcal {P}}_{{{\mathrm {m}}}}^i}(\varvec{u}- \varvec{u}_{{{\mathrm {m}}}}) \big ). \end{aligned}$$

In (28) we denoted by $\varvec{\sigma }^{{\mathcal {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 $\varvec{u}(\xi ,t_{{{\mathrm {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 $\mathop {\bigcup }\limits _i^{N_{{\mathcal {P}}}}{\mathcal {P}}_{{{\mathrm {m}}}}^i$,
2.
[offline] Build at marking time the interpolation operator $\{{\mathcal {J}}_{{\mathcal {P}}_{{{\mathrm {m}}}}^i}\}_{i=1}^{N_{{\mathcal {P}}}}$ from the deformed configuration (by $\varvec{u}(\xi ,t_{{{\mathrm {m}}}})$) to the tag planes $\mathop {\bigcup }\limits _i^{N_{{\mathcal {P}}}}{\mathcal {P}}_{{{\mathrm {m}}}}^i$,
3.
[online] From the estimated displacement, deform the tag planes $\mathop {\bigcup }\limits _i^{N_{{\mathcal {P}}}}{\mathcal {P}}_{{{\mathrm {m}}}}^i$ using (28),
4.
[online] From the estimated (deformed) tag planes $\mathop {\bigcup }\limits _i^{N_{{\mathcal {P}}}}{\widehat{{\mathcal {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 $\varvec{u}(\xi ,t_{{{\mathrm {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

$$\begin{aligned} \begin{pmatrix} 0 &{} K\\ -K&{} -C\end{pmatrix} Y= \lambda \begin{pmatrix} K&{} 0 \\ 0 &{} M\end{pmatrix} Y, \end{aligned}$$

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

$$\begin{aligned} \begin{pmatrix} \gamma K(\varepsilon K+ H^{\intercal }M_{{{\text {obs}}}}H)^{-1}H^\intercal M_{{{\text {obs}}}}H&{} K\\ -K&{} -C\end{pmatrix} Y= \lambda \begin{pmatrix} K&{} 0 \\ 0 &{} M\end{pmatrix} Y, \end{aligned}$$

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 spectral problem reads

$$\begin{aligned} \begin{pmatrix} \gamma K(\varepsilon K+ D_uJ^{\intercal }M_{{{\text {obs}}}}D_uJ)^{-1}D_uJ^\intercal M_{{{\text {obs}}}}D_uJ&{} K\\ -K&{} -C\end{pmatrix} Y= \lambda \begin{pmatrix} K&{} 0 \\ 0 &{} M\end{pmatrix} Y, \end{aligned}$$

where the discrepancy operator corresponds to the one detailed in (13).

Using the optimal criterion on the gain $\gamma $ 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\,{{\mathrm {mm}}}$, to $3.5\,{{\mathrm {mm}}}$ and $0.25\,{{\mathrm {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\,{{\mathrm {mm}}}$. Denoting by $M_{{\mathcal {T}}}$ the surface mass matrix computed on the set of tags, we define the measurement observation norm $S$ by setting, in (23)

$$\begin{aligned} M_{{{\text {obs}}}}= \frac{1}{m} M_{{\mathcal {T}}}. \end{aligned}$$

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 steps— set in our simulations to $2.5\cdot 10^{-4}~$s—and no noise is added. Therefore, following [7] we rescale

$$\begin{aligned} m = \frac{\Delta {{\mathrm {t}}}_{{{\mathrm {obs}}}}}{\Delta {{\mathrm {t}}}} m_{{{\mathrm {obs}}}} = 25~m_{{{\mathrm {obs}}}}, \end{aligned}$$

and set a high confidence in the observations with $m_{{{\mathrm {obs}}}} = (0.65\,{{\mathrm {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 three-dimensional 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–12) and apical (regions 14–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 $\Gamma _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 different—coupled—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.

References

Nash MP, Hunter PJ. Computational mechanics of the heart. J Elast. 2000;61(1–3):113–41.
Article MATH Google Scholar
Sainte-Marie J, Chapelle D, Cimrman R, Sorine M. Modeling and estimation of the cardiac electromechanical activity. Comput Struct. 2006;84(28):1743–59.
Article MathSciNet Google Scholar
Chabiniok R, Wang VY, Hadjicharalambous M, Asner L, Lee J, Sermesant M, et al. Multiphysics and multiscale modelling, data-model fusion and integration of organ physiology in the clinic: ventricular cardiac mechanics. Interface Focus. 2016;6(2):20150083.
Article Google Scholar
Asch M, Bocquet M, Nodet M. Data assimilation: methods, algorithms, and applications. Fundamentals of algorithms. Philadelphia: SIAM; 2016.
Book MATH Google Scholar
Sermesant M, Moireau P, Camara O, Sainte-Marie J, Andriantsimiavona R, Cimrman R, et al. Cardiac function estimation from MRI using a heart model and data assimilation: advances and difficulties. Med Image Anal. 2006;10(4):642–56.
Article Google Scholar
Moireau P, Chapelle D, LeTallec P. Filtering for distributed mechanical systems using position measurements: perspectives in medical imaging. Inverse Probl. 2009;25(3):035010.
Article MathSciNet MATH Google Scholar
Chabiniok R, Moireau P, Lesault PF, Rahmouni A, Deux JF, Chapelle D. Estimation of tissue contractility from cardiac cine-MRI using a biomechanical heart model. Biomech Model Mechanobiol. 2012;11(5):609–30.
Article Google Scholar
Zerhouni EA, Parish D, Rogers W, Yang A, Shapiro E. Human heart: tagging with MR imaging—a method for noninvasive assessment of myocardial motion. Radiology. 1988;169(1):59–63.
Article Google Scholar
Axel L, Dougherty L. MR imaging of motion with spatial modulation of magnetization. Radiology. 1989;171(3):841–5.
Article Google Scholar
Axel L, Montillo A, Kim D. Tagged magnetic resonance imaging of the heart: a survey. Med Image Anal. 2005;9(4):376–93.
Article Google Scholar
Ryf S, Spiegel MA, Gerber M, Boesiger P. Myocardial tagging with 3D CSPAMM. Magn Reson Med. 2002;16(3):320–5.
Google Scholar
Rutz AK, Ryf S, Plein S, Boesiger P, Kozerke S. Accelerated whole-heart 3D CSPAMM for myocardial motion quantification. Magn Reson Med. 2008;59(4):755–63.
Article Google Scholar
Kalman RE, Bucy RS. New results in linear filtering and prediction theory. J Basic Eng. 1961;83(3):95–108.
Article MathSciNet Google Scholar
Bensoussan A. Filtrage Optimal des Systèmes Linéaires. Paris: Dunod; 1971.
MATH Google Scholar
Simon D. Optimal state estimation: Kalman, $H_{\infty }$ and nonlinear approaches. Hoboken: Wiley-Interscience; 2006.
Book Google Scholar
Luenberger DG. An introduction to observers. IEEE Trans Autom Control. 1971;16(6):596–602.
Article Google Scholar
Lakshmivarahan S, Lewis JM. Nudging methods: a critical overview. In: Park SK, Xu L, editors. Data assimilation for atmospheric, oceanic, and hydrologic Applications, vol. XVIII. Berlin: Springer; 2008.
Google Scholar
Moireau P, Chapelle D, Le Tallec P. Joint state and parameter estimation for distributed mechanical systems. Comput Method Appl Mech Eng. 2008;197(6):659–77.
Article MathSciNet MATH Google Scholar
Zhang Q, Clavel A. Adaptive observer with exponential forgetting factor for linear time varying systems. In: Proc. of 40th IEEE conference on decision and control, vol. 4. IEEE; 2001. p. 3886–91.
Pham DT. Stochastic methods for sequential data assimilation in strongly nonlinear systems. J Mar Syst. 2001;129(5):1194–207.
Google Scholar
Moireau P, Bertoglio C, Xiao N, Figueroa CA, Taylor CA, Chapelle D, et al. Sequential identification of boundary support parameters in a fluid-structure vascular model using patient image data. Biomech Model Mechanobiol. 2013;12(3):475–96.
Article Google Scholar
Dokos S, Smaill BH, Young AA, LeGrice IJ. Shear properties of passive ventricular myocardium. Am J Physiol Circ Physiol. 2002;283(6):2650–9.
Article Google Scholar
Guccione JM, Costa KD, McCulloch AD. Finite element stress analysis of left ventricular mechanics in the beating dog heart. J Biomech. 1995;28(10):1167–77.
Article Google Scholar
Vetter FJ, McCulloch AD. Three-dimensional stress and strain in passive rabbit left ventricle: a model study. Ann Biomed Eng. 2000;28(7):781–92.
Article Google Scholar
Caruel M, Chabiniok R, Moireau P, Lecarpentier Y, Chapelle D. Dimensional reductions of a cardiac model for effective validation and calibration. Biomech Model Mechanobiol. 2014;13(4):897–914.
Article Google Scholar
Holzapfel GA, Ogden RW. Constitutive modelling of passive myocardium: a structurally based framework for material characterization. Philos Trans R Soc. 1902;2009(367):3445–75.
MATH Google Scholar
Chapelle D, Le Tallec P, Moireau P, Sorine M. Energy-preserving muscle tissue model: formulation and compatible discretizations. Int J Multiscale Comput Eng. 2012;10(2):189–211.
Article Google Scholar
Bestel J, Clément F, Sorine M. A biomechanical model of muscle contraction. In: Proc. of MICCAI conference. 2001. p. 1159–61.
Guttman M, Prince JL, McVeigh ER. Tag and contour detection in tagged MR images of the left ventricle. IEEE Trans Med Imaging. 1994;13(1):74–88.
Article Google Scholar
Mosher T, Smith M. A DANTE tagging sequence for the evaluation of translational simple motion. Magn Reson Med. 1990;15(2):334–9.
Article Google Scholar
Fischer SE, McKinnon GC, Maier SE, Boesiger P. Improved myocardial tagging contrast. Magn Reson Med. 1993;30(2):191–200.
Article Google Scholar
Shi W, Zhuang X, Wang H, Duckett S, Luong DV, Tobon-Gomez C, et al. A comprehensive cardiac motion estimation framework using both untagged and 3-D tagged MR images based on nonrigid registration. IEEE Trans Med Imaging. 2012;31(6):1263–75.
Article Google Scholar
Kuijerm JPA, Marcus JT, Götte JW, Van Rossum AC, Heethaar RM. Three-dimensional myocardial strain analysis based on short- and long- axis magnetic resonance tagged images using a 1D displacement field. J Magn Reson Imaging. 2000;18(5):553–64.
Article Google Scholar
Chandrashekara R, Mohiaddin RH, Rueckert D. Analysis of 3-D myocardial motion in tagged MR images using nonrigid image registration. IEEE Trans Med Imaging. 2004;23(10):1245–50.
Article Google Scholar
Declerck J, Ayache N, McVeigh ER. Use of a 4D planispheric transformation for the tracking and analysis of LV motion with tagged MR images. In: Proc. of SPIE conference on medical imaging. 1999. p. 69–80.
Kerwin WS, Prince JL. MR tag surface tracking using a spatio-temporal filter/interpolator. In: In Proc. ICIP-98, vol. 1. 1998. p. 699–703.
Amini AA, Chen Y, Elayyadi M. Tag surface reconstruction and tracking of myocardial beads from SPAMM-MRI with parametric B-spline surfaces. IEEE Trans Med Imaging. 2001;20(2):94–103.
Article Google Scholar
Amini AA, Chen Y, Curwen RW, Mani V, Sun J. Coupled B-Snake grids and constrained thin-plate splines for analysis of 2D tissue deformations from tagged MRI. IEEE Trans Med Imaging. 1998;17(3):344–56.
Article Google Scholar
Osman NF, Prince JL. Direct calculation of 2D components of myocardial strain using sinusoidal MR tagging. In: Proc. SPIE Med. Imag. Conf. 1998. p. 142–52.
Clarysse P, Basset C, Khouas L, Croisille P, Friboulet D, Odet C, et al. Two-dimensional spatial and temporal displacement and deformation field fitting from cardiac magnetic resonance tagging. Med Image Anal. 2000;4(3):253–68.
Article Google Scholar
Ledesma-Carbayo MJ, Bajo A, Santa Marta C, Perez-David E, Garcia-Fernandez MA, Desco M, et al. Fully automatic cardiac motion estimation for tagged MRI using non-rigid registration techniques. In: Proc. of computers in cardiology. IEEE. 2006. p. 305–8.
Bruurmijn LCM, Kause HB, Filatova OG, Duits R, Fuster A, Florack LMJ, et al. Myocardial deformation from local frequency estimation in tagging MRI. In: Proc. of FIMH conference. Springer; 2013. p. 284–91.
Stoeck CT, Manka R, Boesiger P, Kozerke S. Undersampled cine 3D tagging for rapid assessment of cardiac motion. J Cardiovasc Magn Reson. 2012;14(1):60.
Article Google Scholar
O’Dell W, Moore CC, Hunter WC, Zerhouni EA, McVeigh ER. Three-dimensional myocardial deformations: calculation with displacement field fitting to tagged MR Images. Radiology. 1995;195(3):829–35.
Article Google Scholar
Pan L, Prince JL, Lima ACJ, Osman NF. Fast tracking of cardiac motion using 3D-HARP. IEEE Trans Biomed Eng. 2005;52(8):1425–35.
Article Google Scholar
Amini AA, Shi P, Constable RT. Energy-minimizing deformable grids for tracking tagged MR cardiac images. In: Computers in cardiology; 1992. p. 651–4.
Radeva P, Amini AA, Huang J. Deformable B-Solids and implicit snakes for 3D localization and tracking of SPAMM MRI-Data. In: Proc. of the Workshop on mathematical methods in biomedical image analysis. IEEE. 1996. p. 192–201.
Rosset A, Spadola L, Ratib O. OsiriX: an open-source software for navigating in multidimensional DICOM images. J Digit Imaging. 2004;17(3):205–16.
Article Google Scholar
Prince JL, McVeigh ER. Motion estimation from tagged MR image sequences. IEEE Trans Med Imaging. 1992;11(2):238–49.
Article Google Scholar
Dougherty L, Asmuth JC, Blom AS, Axel L, Kumar R. Validation of an optical flow method for tag displacement estimation. IEEE Trans Med Imaging. 1999;18(4):359–63.
Article Google Scholar
Osman NF, Kerwin WS, McVeigh ER, Prince JL. Cardiac motion tracking using CINE harmonic phase (HARP) magnetic resonance imaging. Magn Reson Med. 1999;42(6):1048–60.
Article Google Scholar
Osman NF, Prince JL. Visualizing myocardial function using HARP MRI. Phys Med Biol. 2000;45(6):1665–82.
Article Google Scholar
Garcia-Barnés J, Gil D, Pujadas S, Carreras F. A variational framework for assessment of the left ventricle motion. Math Model Nat Phenom. 2008;3(06):76–100.
Article MathSciNet MATH Google Scholar
Kalman R. A new approach to linear filtering and prediction problems. Trans ASME J Basic Eng. 1960;82:35–45.
Article MathSciNet Google Scholar
Moireau P, Chapelle D. Reduced-order unscented Kalman filtering with application to parameter identification in large-dimensional systems. ESAIM COCV. 2011;17(02):380–405.
Article MathSciNet MATH Google Scholar
Xi J, Lamata P, Lee J, Moireau P, Chapelle D, Smith N. Myocardial transversely isotropic material parameter estimation from in-silico measurements based on a reduced-order unscented Kalman filter. J Mechan Behav Biomed Mater. 2011;4(7):1090–102.
Article Google Scholar
Marchesseau S, Delingette H, Sermesant M, Cabrera-Lozoya R, Tobon-Gomez C, Moireau P, et al. Personalization of a cardiac electromechanical model using reduced order unscented Kalman filtering from regional volumes. Med Image Anal. 2013;17(7):816–29.
Article Google Scholar
Julier SJ. Reduced sigma point filters for the propagation of means and covariances through nonlinear transformations. In: Proc. of the American control conference, vol. 2. 2002. p. 887–92.
Cerqueira MD, Weissman NJ, Dilsizian V, Jacobs AK, Kaul S, Laskey WK, et al. Standardized myocardial segmentation and nomenclature for tomographic imaging of the heart: A statement for healthcare professionals from the cardiac imaging committee of the council on clinical cardiology of the American Heart Association. Circulation. 2002;105(4):539–42.
Article Google Scholar
Chen C, Qin C, Qiu H, Tarroni G, Duan J, Bai W, et al. Deep learning for cardiac image segmentation: a review. Front Cardiovasc Med. 2020;7:25.
Article Google Scholar
Toussaint N, Stoeck CT, Schaeffter T, Kozerke S, Sermesant M, Batchelor PG. In vivo human cardiac fibre architecture estimation using shape-based diffusion tensor processing. Med Image Anal. 2013;17(8):1243–55.
Article Google Scholar
Doste R, Soto-Iglesias D, Bernardino G, Alcaine A, Sebastian R, Giffard-Roisin S, et al. A rule-based method to model myocardial fiber orientation in cardiac biventricular geometries with outflow tracts. Int J Numer Methods Biomed Eng. 2019;35(4):e3185 E3185 cnm.3185.
Article Google Scholar
Nagler A, Bertoglio C, Stoeck CT, Kozerke S, Wall WA. Maximum likelihood estimation of cardiac fiber bundle orientation from arbitrarily spaced diffusion weighted images. Med Image Anal. 2017;39:56–77.
Article Google Scholar
Rausch MK, Genet M, Humphrey JD. An augmented iterative method for identifying a stress-free reference configuration in image-based biomechanical modeling. J Biomech. 2017;58:227–31.
Article Google Scholar
Miyashita H. Clinical assessment of central blood Pressure. Curr Hypertens Rev. 2012;8(2):80–90.
Article Google Scholar
Laleg TM, Crépeau E, Sorine M. Separation of arterial pressure into a nonlinear superposition of solitary waves and a windkessel flow. Biomed Signal Process Control. 2007;2(3):163–70.
Article Google Scholar
Bertoglio C, Barber D, Gaddum N, Valverde I, Rutten M, Beerbaum P, et al. Identification of artery wall stiffness: in vitro validation and in vivo results of a data assimilation procedure applied to a 3D fluid-structure interaction model. J Biomech. 2014;47(5):1027–34.
Article Google Scholar
Chapelle D, Fragu M, Mallet V, Moireau P. Fundamental principles of data assimilation underlying the Verdandi library: applications to biophysical model personalization within euHeart. Med Biol Eng Comput. 2013;51(11):1221–33.
Article Google Scholar
Vaillant M, Glaunès J. Surface matching via currents. In: Golland P, Fischl B, editors. Information processing in medical imaging. Berlin: Springer; 2005. p. 381–92.
Chapter Google Scholar
Durrleman S, Pennec X, Trouvé A, Thompson P, Ayache N. Inferring brain variability from diffeomorphic deformations of currents: an integrative approach. Med Image Anal. 2008;12(5):626–37.
Article Google Scholar
Younes L. Shapes and diffeomorphisms Applied mathematical sciences. Berlin: Springer; 2010.
Book MATH Google Scholar
Dautray R, Lions JL. Mathematical analysis and numerical methods for science and technology. Evolution problems I, vol. 5. Berlin: Springer; 1992.
MATH Google Scholar
Chapelle D, Cîndea N, de Buhan M, Moireau P. Exponential convergence of an observer based on partial field measurements for the wave equation. Math Prob Eng. 2012;2012:1–12.
Article MathSciNet MATH Google Scholar
Liu K. Locally distributed control and damping for the conservative systems. SIAM J Control Optim. 1997;35(5):1574–90.
Article MathSciNet MATH Google Scholar
Daoulatli M, Dehman B, Khenissi M. Local energy decay for the elastic system with nonlinear damping in an exterior Domain. SIAM J Control Optim. 2010;48(8):5254–75.
Article MathSciNet MATH Google Scholar
Bardos C, Lebeau G, Rauch J. Un exemple d’utilisation des notions de propagation pour le contrôle et la stabilisation des problèmes hyperboliques. Rendiconti del Seminario Matematico del Universita Politecnico Torino. Fascicolo speciale(Hyperbolic Equations (1987)). 1988. pp. 12–31.

Download references

Funding

The authors gratefully acknowledge support from the European Union’s Seventh Framework Programme for research, technological development and demonstration, under Grant agreement #611823 (VP2HF Project).

Author information

Authors and Affiliations

CEA List, CEA, 91190, Saclay, France
Alexandre Imperiale
Project-Team M3DISIM, Inria Saclay-Ile-de-France, Inria, 91128, Palaiseau, France
Alexandre Imperiale, Dominique Chapelle & Philippe Moireau
LMS, Ecole Polytechnique, CNRS, Institut Polytechnique de Paris, 91128, Palaiseau, France
Dominique Chapelle & Philippe Moireau

Authors

Alexandre Imperiale
View author publications
You can also search for this author in PubMed Google Scholar
Dominique Chapelle
View author publications
You can also search for this author in PubMed Google Scholar
Philippe Moireau
View author publications
You can also search for this author in PubMed Google Scholar

Contributions

AI, DC and PM contributed to the conception, design, analysis, interpretations’, and the drafting and revising of the manuscript. Moreover, AI and PM contributed to the software engineering of the method. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Philippe Moireau.

Ethics declarations

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendices

Appendix A: The illuminating example of efficient state feedback in elastodynamics

In the following very simplified example, we aim at mathematically proving why a simple state feedback can have a comparable efficiency with respect to a Kalman-based feedback in the context of elastodynamics problems. Let us then consider a simplified configuration where the model and the observation operator are linear. More precisely, we assume a linear elastic system in which the Cauchy stress tensor is given by

$$\begin{aligned} \varvec{\sigma }= \varvec{A}: \varvec{\varepsilon }(\varvec{u}), \end{aligned}$$

where the elasticity tensor $\varvec{A}$ is assumed to be constant and isotropic. Defining $\Delta _e(\cdot ) = \mathop {\mathbf {div}}( \varvec{A}: \varvec{\varepsilon }(\cdot ) )$, our model simply reads

$$\begin{aligned}&\left\{ \begin{array}{ll} {\dot{\varvec{u}}} = \varvec{v}, &{} \quad \text{ in } \Omega _0\times (0,T) \\ \rho {\dot{\varvec{v}}} - \Delta _e \varvec{u}= \varvec{f},&{} \quad \text{ in } \Omega _0\times (0,T) \\ \varvec{u}= 0, &{} \quad {{\text {on }}} \partial \Omega _0\times (0,T). \end{array} \right. \end{aligned}$$

(29)

The external load is a time-dependent regular function $\varvec{f} \in C^1([0,T],{\mathcal {L}}^2(\Omega _0))$. We introduce ${\mathcal {Y}}^{v}= {\mathcal {L}}^2(\Omega _0)^3$, the displacement space ${\mathcal {Y}}^{u}= {\mathcal {H}}^1_0(\Omega _0)^3$, and $ {\mathcal {Y}}= {\mathcal {Y}}^{u}\times {\mathcal {Y}}^{v}$. Using the Korn and Poincaré inequalities, ${\mathcal {Y}}^{u}$ is an Hilbert space with the following scalar product

$$\begin{aligned} \forall (\varvec{u}_1,\varvec{u}_2) \in {\mathcal {Y}}^{u}, \quad (\varvec{u}_1, \varvec{u}_2 )_{{\mathcal {E}}_0} = \int _{\Omega _0} \varvec{\varepsilon }(\varvec{u}_1) : \varvec{A}: \varvec{\varepsilon }(\varvec{u}_2) \, {{\mathrm {d}}}\Omega . \end{aligned}$$

In a semi-group theory context we introduce the semi-group generator ${{\mathrm {A}}}\in {\mathcal {L}}({\mathcal {D}}({{\mathrm {A}}}),{\mathcal {Y}})$ with

$$\begin{aligned} {{\mathrm {A}}}= \begin{pmatrix} 0 &{} \mathbb {1} \\ \tfrac{1}{\rho } \Delta _e &{} 0 \end{pmatrix}, \end{aligned}$$

and we can prove that (29) admits a classical solution $C^0([0,T],{\mathcal {Y}}^{u}) \cap C^1([0,T],{\mathcal {Y}}^{v})$ for every initial condition in the domain [72]

$$\begin{aligned} {\mathcal {D}}({{\mathrm {A}}}) = \left\{ (\varvec{u},\varvec{v}) \in {\mathcal {Y}}^{u}\times {\mathcal {Y}}^{v}, \quad \mathop {\mathbf {div}}( \varvec{A}: \varvec{\varepsilon }(\varvec{u}) ) \in {\mathcal {Y}}^{v}\right\} . \end{aligned}$$

Moreover, the operator ${{\mathrm {A}}}$ is skew-adjoint, implying that

$$\begin{aligned} \forall t > 0, \quad \left\| {{\mathrm {y}}}(t)\right\| _{{\mathcal {Y}}} = \left\| {{\mathrm {y}}}(0)\right\| _{{\mathcal {Y}}}, \end{aligned}$$

corresponding to the energy balance on the system (29) with ${\mathcal {E}} = \frac{1}{2} \left\| {{\mathrm {y}}}\right\| _{{\mathcal {Y}}}^2$. Using the semi-group theory, we rewrite the dynamics of the model in the abstract state-space form

$$\begin{aligned} {\dot{{{\mathrm {y}}}}} = {{\mathrm {A}}}{{\mathrm {y}}}+ {{\mathrm {R}}}, \quad {{\mathrm {y}}}(0)= {{\mathrm {y}}}_0+ \zeta _{{{\mathrm {y}}}}. \end{aligned}$$

Concerning this model, we assume that we have at our disposal some measurements of the displacements. We introduce the observation operator ${{\mathrm {H}}}= ({{\mathrm {H}}}_{|{\mathcal {Y}}^{u}} \quad 0)$ such that

where $\omega _0$ is an internal open subdomain of $\Omega _0$ and ${\mathcal {Z}}= {\mathcal {H}}^1(\omega _0)^3$. Here ${{\mathrm {H}}}$ and ${{\mathrm {H}}}_{|{\mathcal {Y}}^{u}}$ as ${{\mathrm {H}}}= ({{\mathrm {H}}}_{|{\mathcal {Y}}^{u}} \quad 0)$ correspond to the same operators with different input spaces.

Using the extension $\varvec{\psi } = \mathbf {Ext}_{\omega _0}(\varvec{\varphi })$ defined by

$$\begin{aligned} {\left\{ \begin{array}{ll} - \Delta _e(\varvec{\psi }) = 0, &{} \quad {{\text {in}}} ~\Omega _0\\ \varvec{\psi } = \varvec{\varphi }, &{}\quad {{\text {in}}} ~\omega _0\\ \varvec{\psi } = 0, &{}\quad {{\text {on}}} ~\partial \Omega _0, \end{array}\right. } \end{aligned}$$

(30)

we can prove the following property.

Proposition A.1

For any $(\varvec{\varphi }_1,\varvec{\varphi }_2) \in {\mathcal {Z}}^2$, the bilinear form $(\mathbf {Ext}_{\omega _0}(\varvec{\varphi }_1),\mathbf {Ext}_{\omega _0}(\varvec{\varphi }_2))_{{\mathcal {E}}_0}$ defines a scalar product on ${\mathcal {Z}}= {\mathcal {H}}^1(\omega _0)^3$.

Proof

The proof is a simple extension of the property proven by [73] for scalar equations. Let $\varvec{\varphi }$ be an element of ${\mathcal {Z}}$. The only difficulty lies in proving the norm equivalence with $\left\| \varvec{\varphi }\right\| _{{\mathcal {H}}^1(\omega _0)^3}^2$. First, we have

$$\begin{aligned} \left\| \varvec{\varphi }\right\| _{{\mathcal {H}}^1(\omega _0)^3}^2&= \left\| \varvec{\nabla }\,\varvec{\varphi }\right\| _{{\mathcal {L}}^2(\omega _0)^3}^2 + \left\| \varvec{\varphi }\right\| _{{\mathcal {L}}^2(\omega _0)^3}^2 \\&\le \left\| \varvec{\nabla }\, \mathbf {Ext}_{\omega _0}(\varvec{\varphi })\right\| _{{\mathcal {L}}^2(\Omega _0)^3}^2 +\left\| \mathbf {Ext}_{\Omega _0}(\varvec{\varphi })\right\| _{{\mathcal {L}}^2(\Omega _0)^3}^2\\&\le (1+C_p) \left\| \varvec{\nabla }\, \mathbf {Ext}_{\omega _0}(\varvec{\varphi })\right\| _{{\mathcal {L}}^2(\Omega _0)^3}^2\\&\le C_k(1+C_p) \left\| \mathbf {Ext}_{\omega _0}(\varvec{\varphi })\right\| _{{\mathcal {E}}_0}^2, \end{aligned}$$

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 $\Omega _0\backslash \omega _0$ with respect to the data, there exists a constant $C_d >0$ such that for any $\varvec{\psi } = \mathbf {Ext}_{\omega _0}(\varvec{\varphi }) $ we have

$$\begin{aligned} \int _{\Omega _0\backslash \omega _0} \varvec{\varepsilon }(\varvec{\psi }) : \varvec{A}: \varvec{\varepsilon }(\varvec{\psi }) \, {{\mathrm {d}}}\Omega \le C_d \Vert \varvec{\varphi }_{|\partial \omega _0}\Vert _{{\mathcal {H}}^{\frac{1}{2}}(\partial \omega _0)^3}^2. \end{aligned}$$

Hence, denoting by $C_t$ the constant arising from the continuity of the trace operator, we have

$$\begin{aligned} \left\| \mathbf {Ext}_{\omega _0}(\varvec{\varphi })\right\| _{{\mathcal {E}}_0}^2&\le \int _{\omega _0} \varvec{\varepsilon }(\varvec{\psi }) : \varvec{A}: \varvec{\varepsilon }(\varvec{\psi }) \, {{\mathrm {d}}}\Omega + C_d \left\| \varvec{\varphi }\right\| _{{\mathcal {H}}^{\frac{1}{2}}(\partial \omega _0)^3}^2 \\&\le C_a \left\| \varvec{\nabla }\, \varvec{\varphi }\right\| _{{\mathcal {L}}^{2}(\omega _0)^3}^2 + C_d \Vert \varvec{\varphi }_{|\partial \omega _0}\Vert ^2_{{\mathcal {H}}^{\frac{1}{2}}(\partial \omega _0)^3} \\&\le (C_a + C_d C_t) \Vert \varvec{\varphi } \Vert ^2_{{\mathcal {H}}^{1}( \omega _0)^3}, \end{aligned}$$

which completes the proof. $\square $

It is now possible to define the adjoint of the observation operator.

Proposition A.2

The operator ${{\mathrm {H}}}$ is bounded from ${\mathcal {Y}}$ to ${\mathcal {Z}}$ and ${{\mathrm {H}}}^*$ is given by

$$\begin{aligned} {{\mathrm {H}}}^* = \begin{pmatrix} {{\mathrm {H}}}_{|{\mathcal {Y}}^{u}}^* \\ 0 \end{pmatrix} {{\text { with }}} {{\mathrm {H}}}_{|{\mathcal {Y}}^{u}}^*\, :\, \left| \begin{aligned} {\mathcal {Z}}&\rightarrow {\mathcal {Y}}^{u}\\ \varvec{\varphi }&\mapsto \mathbf {Ext}_{\omega _0}(\varvec{\varphi }) \end{aligned} \right. \end{aligned}$$

Proof

Let us first prove that ${{\mathrm {H}}}$ is bounded. We consider $\varvec{\psi }\in {\mathcal {Y}}^{u}$ and $\varvec{\varphi }$ such that $\varvec{\varphi } = {{\mathrm {H}}}_{|{\mathcal {Y}}^{u}}\,\varvec{\psi }$. We have directly, from norm equivalences,

$$\begin{aligned} \left\| \varvec{\varphi }\right\| _{\mathcal {Z}}^2 = \left\| \mathbf {Ext}_{\omega _0}(\varvec{\varphi })\right\| _{{\mathcal {E}}_0}^2 \le C_1 \left\| \varvec{\varphi }\right\| _{{\mathcal {H}}^1(\omega _0)^3}^2 \le C_1 \left\| \varvec{\psi }\right\| _{{\mathcal {H}}^1_0(\Omega _0)^3}^2. \end{aligned}$$

Then, we have that for all $\varvec{\varphi } \in {\mathcal {Z}}$ and $\varvec{v}^{\sharp }\in {\mathcal {Y}}^{u}$

$$\begin{aligned} (\varvec{\varphi }, {{\mathrm {H}}}_{|{\mathcal {Y}}^{u}}\, \varvec{v}^{\sharp })_{{\mathcal {Z}}} = \int _{\Omega _0} \varvec{\varepsilon }(\mathbf {Ext}_{\omega _0}(\varvec{\varphi })) : \varvec{A}: \varvec{\varepsilon }(\mathbf {Ext}_{\omega _0}(\varvec{v}^{\sharp }_{|\omega _0})) \, {{\mathrm {d}}}\Omega . \end{aligned}$$

By the variational characterization of the extension (5) we have

$$\begin{aligned} \int _{\Omega _0} \varvec{\varepsilon }(\mathbf {Ext}_{\omega _0}(\varvec{\varphi })) : \varvec{A}: \varvec{\varepsilon }(\mathbf {Ext}_{\omega _0}(\varvec{v}^{\sharp }_{|\omega _0})) \, {{\mathrm {d}}}\Omega = \int _{\Omega _0} \varvec{\varepsilon }(\mathbf {Ext}_{\omega _0}(\varvec{\varphi })) : \varvec{A}: \varvec{\varepsilon }(\varvec{v}^{\sharp }) \, {{\mathrm {d}}}\Omega , \end{aligned}$$

since $\varvec{v}^{\sharp }_{|\omega _0} - \mathbf {Ext}_{\omega _0}(\varvec{v}^{\sharp }_{|\omega _0}) = 0$ on $\omega _0$. Therefore $(\varvec{\varphi }, {{\mathrm {H}}}_{|{\mathcal {Y}}^{u}}\, \varvec{v}^{\sharp })_{{\mathcal {Z}}} = ({{\mathrm {H}}}_{|{\mathcal {Y}}^{u}}^* \varvec{\varphi }, \varvec{v}^{\sharp })_{{\mathcal {E}}_0}, $ and ${{\mathrm {H}}}^*$ is given by

$$\begin{aligned} {{\mathrm {H}}}^* : \left| \begin{aligned} {\mathcal {Z}}&\rightarrow {\mathcal {Y}}^{u}\\ \varvec{\varphi }&\mapsto \begin{pmatrix}{{\text {Ext}}}_{\omega _0}(\varvec{\varphi })\\ 0\end{pmatrix}. \end{aligned} \right. \end{aligned}$$

$\square $

We can now define the observer by the dynamics

which in strong form reads

(31)

which converges to the solution of (29) under the observability condition given by the next theorem, see e.g. [74] for a proof.

Theorem A.3

Let ${{\mathrm {A}}}$ be a time-independent skew-adjoint operator generating a group and ${{\mathrm {H}}}\in {\mathcal {L}}({\mathcal {Y}},{\mathcal {Z}})$. The error system ${\widetilde{{{\mathrm {y}}}}}$ of dynamics

$$\begin{aligned} \dot{{\widetilde{{{\mathrm {y}}}}}} = ({{\mathrm {A}}}- \gamma {{\mathrm {H}}}^* {{\mathrm {H}}}) {\widetilde{{{\mathrm {y}}}}} \end{aligned}$$

is exponentially stable if the following observability condition is satisfied: there exists two constants $(C_{{{\text {st}}}},T)$ such that for any solution ${{\mathrm {y}}}$ of ${\dot{{{\mathrm {y}}}}} = {{\mathrm {A}}}{{\mathrm {y}}}$, we have

$$\begin{aligned} \int _0^T \left\| {{\mathrm {H}}}{{\mathrm {y}}}(s)\right\| ^2_{\mathcal {Z}}\, {{\mathrm {d}}}s \ge C_{{{\text {st}}}} \left\| {{\mathrm {y}}}(0)\right\| ^2_{{\mathcal {Y}}}. \end{aligned}$$

(32)

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_{{\text {st}}}$ and a time T such that every solution of

$$\begin{aligned}&\left\{ \begin{array}{ll} {\dot{\varvec{u}}} = \varvec{v}, &{} \quad \text{ in } \Omega _0\\ \rho {\dot{\varvec{v}}} - \Delta _e \varvec{u}= 0,&{} \quad \text{ in } \Omega _0\\ \varvec{u}= 0, &{} \quad {{\text {on }}} \partial \Omega _0\end{array} \right. \end{aligned}$$

satisfies the observability condition

$$\begin{aligned} \int _{0}^T&\left\| \mathbf {Ext}_{\omega _0}(\mathbb {1}_{\omega _0} \varvec{u})\right\| _{{\mathcal {E}}_0}^2 \, {{\mathrm {d}}}t \ge c^{{\text {st}}} \Bigl ( \left\| \varvec{u}(0)\right\| _{{\mathcal {E}}_0}^2 + \left\| \varvec{v}(0)\right\| _{{\mathcal {L}}^2}^2 \Bigr ), \end{aligned}$$

(33)

then, in the absence of observation error, the observer given by the dynamics (31) converges to the solution ${{\mathrm {y}}}^{{{\text {ref}}}}$ of (29) such that

$$\begin{aligned} {{\mathrm {z}}}= {{\mathrm {H}}}{{\mathrm {y}}}^{{{\text {ref}}}}. \end{aligned}$$

Proof

We have defined the reference trajectory as the solution of

$$\begin{aligned} {\dot{{{\mathrm {y}}}}}^{{{\text {ref}}}} = {{\mathrm {A}}}{{\mathrm {y}}}^{{{\text {ref}}}} + {{\mathrm {R}}}, \end{aligned}$$

and the observer as the solution of

The error ${\widetilde{{{\mathrm {y}}}}} = {{\mathrm {y}}}^{{{\text {ref}}}} - {\widehat{{{\mathrm {y}}}}}$ is then solution of

which, from Theorem A.3 of Appendix A, converges exponentially to 0 for every initial condition when the observability condition (33) is verified. $\square $

Following [75], we define the elastic geometric control condition:

Definition A.1

(Elastic Geometric Control Condition) The elastic geometric control condition is satisfied if every combination of pressure (P) and shear (S) waves ray encounters—in 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 $\omega _0$ as soon as the elastic geometric control condition is satisfied with ${\check{\omega }}_0$ and ${{\text {dist}}}(\Omega _0\backslash \omega _0,{\check{\omega }}_0)>0$.

Proof

For technical reasons, we assume that the elastic geometric control condition is satisfied for an observation domain ${\check{\omega }}_0$ slightly smaller than $\omega _0$, namely, with ${\check{\omega }}_0\subset \omega _0$ and ${{\text {dist}}}(\Omega _0\backslash \omega _0,{\check{\omega }}_0)>0$. We first recall the classical observability result when the velocity is observed. In fact there exists a constant $C_{{{\text {st}}}}$ and a time T such that every solution of (29) satisfies the observability condition

$$\begin{aligned} \int _{0}^{\breve{T}} \left\| \varvec{v}\right\| _{{\mathcal {L}}^2({\check{\omega }}_0)^3}^2 \, {{\mathrm {d}}}t\ge C_{{{\text {st}}}}\Bigl ( \left\| \varvec{u}(0)\right\| _{{\mathcal {E}}_0}^2 + \left\| \varvec{v}(0)\right\| _{{\mathcal {L}}^2(\Omega _0)^3}^2 \Bigr ), \end{aligned}$$

(34)

with $\breve{T}=T-\delta $ for $\delta >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 $\psi \in C_c^\infty ({\overline{\Omega }}_0)$ be a cutoff function satisfying

$$\begin{aligned} \psi (\varvec{\xi }) = \left\{ \begin{array}{l} 0, \quad \text {if } \varvec{\xi }\in \Omega _0\!\backslash \! \omega _0\\ 1, \quad \text {if } \varvec{\xi }\in \breve{\omega }_0 \end{array} \right. \end{aligned}$$

and $0 \le \psi (\varvec{\xi }) \le 1$ for every $\varvec{\xi }\in {\overline{\Omega }}_0$. Denote also $\phi (t) = t^2 ({\breve{T}} - t)^2$. Then, by repeated integrations by parts we obtain

$$\begin{aligned} 0&= \int _0^{\breve{T}} \int _{\omega _0} \phi \psi (\ddot{\varvec{u}} - \Delta _e \varvec{u}) \cdot \varvec{u}\, {{\mathrm {d}}}\Omega \, {{\mathrm {d}}}t \\&= \int _0^{\breve{T}} \int _{\omega _0} \ddot{\phi } \psi \frac{|\varvec{u}|^2}{2} \, {{\mathrm {d}}}\Omega \, {{\mathrm {d}}}t - \int _0^{\breve{T}} \int _{\omega _0} \phi \psi |{\dot{\varvec{u}}}|^2 \, {{\mathrm {d}}}\Omega \, {{\mathrm {d}}}t \\&\quad - \int _0^{\breve{T}} \int _{\omega _0} \phi \, \varvec{\varepsilon }(\varvec{u}) : \varvec{A}: (\varvec{\nabla } \psi \otimes \varvec{u})\, {{\mathrm {d}}}\Omega \, {{\mathrm {d}}}t \\&\quad + \int _0^{\breve{T}} \int _{\omega _0} \phi \psi \, \varvec{\varepsilon }(\varvec{u}) : \varvec{A}: \varvec{\varepsilon }(\varvec{u}) \, {{\mathrm {d}}}\Omega \, {{\mathrm {d}}}t. \end{aligned}$$

Moreover,

$$\begin{aligned}&\int _{\omega _0} \varvec{\varepsilon }(\varvec{u}) : \varvec{A}: (\varvec{\nabla } \psi \otimes \varvec{u}) \, {{\mathrm {d}}}\Omega \\&\quad \le C_{{{\text {st}}}} \left\| \psi \right\| _{{\mathcal {W}}^{1,\infty }} \left\| \varvec{\varepsilon }(\varvec{u})\right\| _{{\mathcal {L}}^2(\omega _0)^3} \left\| \varvec{u}\right\| _{{\mathcal {L}}^2(\omega _0)^3} \\&\quad \le C_{{{\text {st}}}} \left\| \psi \right\| _{{\mathcal {W}}^{1,\infty }} \left\| \varvec{u}\right\| _{{\mathcal {H}}^1(\omega _0)^3}^2. \end{aligned}$$

where $C_{{{\text {st}}}}$ represents a different constant in each line. This identity combined with the properties of the cutoff functions $\phi $ and $\psi $ provides, for any strictly positive $\varepsilon $, the existence of a constant $C_{{{\text {st}}}} > 0$ such that

$$\begin{aligned} \int _\varepsilon ^{\breve{T}-\varepsilon } \int _{{\check{\omega }}_0} |{\dot{\varvec{u}}}|^2 \, {{\mathrm {d}}}\Omega \, {{\mathrm {d}}}t \le C_{{{\text {st}}}} \int _0^{\breve{T}} \int _{\omega _0} \Vert \varvec{u}(\cdot ,t)\Vert _{{\mathcal {H}}^1(\omega _0)^3}^2 \, {{\mathrm {d}}}t . \end{aligned}$$

Substituting ${\breve{T}}+2\varepsilon $ for ${\breve{T}}$ in all the above computations gives

$$\begin{aligned} \int _\varepsilon ^{{\breve{T}} + \varepsilon } \int _{\breve{\omega }_0} |{\dot{\varvec{u}}}|^2\, {{\mathrm {d}}}\Omega \, {{\mathrm {d}}}t \le C_{{{\text {st}}}}\int _0^{{\breve{T}}+2\varepsilon } \Vert \varvec{u}(\cdot ,t)\Vert _{{\mathcal {H}}^1(\omega _0)^3}^2 \, {{\mathrm {d}}}t. \end{aligned}$$

We proceed by making the change of variable $\tau = t - \varepsilon $ in the left-hand side integral, yielding

$$\begin{aligned} \int _0^{\breve{T}} \int _{\breve{\omega }_0} |{\dot{\varvec{u}}}(\varvec{\xi }, \tau + \varepsilon )|^2 \, {{\mathrm {d}}}\Omega \, {{\mathrm {d}}}\tau \le C_{{{\text {st}}}} \int _0^{{\breve{T}}+2\varepsilon } \Vert \varvec{u}\Vert _{{\mathcal {H}}^1(\omega _0)^3}^2 \, {{\mathrm {d}}}t. \end{aligned}$$

(35)

Noting that $\varvec{u}(\varvec{\xi }, t + \varepsilon )$ satisfies the elastodynamics system with initial data $(\varvec{u}(\varvec{\xi },\varepsilon ),{\dot{\varvec{u}}}(\varvec{\xi },\varepsilon ))$ and applying (34) with this shifted solution, we obtain that there exists also $C_{{{\text {st}}}}$ such that

$$\begin{aligned} \int _0^{\breve{T}} \int _{\breve{\omega }_0} |{\dot{\varvec{u}}}(\varvec{\xi }, \tau + \varepsilon )|^2 \, {{\mathrm {d}}}\Omega \, \, {{\mathrm {d}}}\tau \ge C_{{{\text {st}}}} \left( \Vert \varvec{u}(\varepsilon )\Vert ^2_{{\mathcal {E}}_0} + \Vert {\dot{\varvec{u}}}(\varepsilon )\Vert ^2_{{\mathcal {L}}^2(\Omega _0)^3} \right) . \end{aligned}$$

(36)

Combining (35), (36) and the fact that the energy of the solution of the elastodynamics equation is exactly conserved over time, we have our observability inequality (33) upon choosing $\varepsilon =\frac{\delta }{2}$. $\square $

Appendix B: Analysis of the prediction–correction scheme

This appendix is dedicated to the analysis of the prediction–correction scheme proposed in “Time discretization of the state observer” section. More specifically, we are interested in showing that the linearization procedure applied in the context of nonlinear observation operators induces a dissipative behavior for the linearized time-discrete estimation error. This property is a crucial aspect when building relevant state observers.

Linear model and linear observation operator

To start with, we assume that the dynamical system satisfied by the target trajectory is given by

$$\begin{aligned} {\left\{ \begin{array}{ll} {\dot{{{\mathrm {y}}}}}(t) = ({{\mathrm {A}}}+ \upeta {{\mathrm {V}}}){{\mathrm {y}}}(t)\\ {{\mathrm {y}}}(0) = {{\mathrm {y}}}_0 + \zeta _{{\mathrm {y}}}, \end{array}\right. } \end{aligned}$$

(37)

where ${{\mathrm {A}}}$ is a skew-adjoint operator, ${{\mathrm {V}}}$ is a self-adjoint and semi-negative operator and $\upeta \ge 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 $\Delta {{\mathrm {t}}}$ the (constant) time step of the numerical procedure, the time-discrete observations read ${{\mathrm {z}}}^n = {{\mathrm {H}}}{{\mathrm {y}}}(n\Delta {{\mathrm {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

$$\begin{aligned} {\widetilde{{{\mathrm {y}}}}}^{n}_+ = {{\mathrm {y}}}(n\Delta {{\mathrm {t}}}) - {\widehat{{{\mathrm {y}}}}}^{n}_+, \end{aligned}$$

(39)

and associating a prediction error with

$$\begin{aligned} {\widetilde{{{\mathrm {y}}}}}^{n+1}_- = {{\mathrm {y}}}((n+1)\Delta {{\mathrm {t}}}) - {\widehat{{{\mathrm {y}}}}}^{n+1}_-, \end{aligned}$$

(40)

we can determine the time-discrete dynamics satisfied by the estimation error.

Proposition B.1

Assuming that ${{\mathrm {y}}}\in C^3([0,T],{\mathcal {Y}})$, then the estimation error satisfies the following discrete dynamical system

with $\varepsilon ^n=O\bigl (\Delta {{\mathrm {t}}}^2\left\| \dddot{{{\mathrm {y}}}}\right\| _{C^3([0,T],{\mathcal {Y}})}\bigr )$.

Proof

(41b) is directly inferred from the definition of the prediction estimation error and using (38b), namely,

$$\begin{aligned} \begin{aligned} {\widetilde{{{\mathrm {y}}}}}^{n+1}_-&= {{\mathrm {y}}}((n+1)\Delta {{\mathrm {t}}}) - {\widehat{{{\mathrm {y}}}}}^{n+1}_- \\&= {{\mathrm {y}}}((n+1)\Delta {{\mathrm {t}}}) - {\widehat{{{\mathrm {y}}}}}^{n+1}_+ + \gamma \Delta {{\mathrm {t}}}{{{\mathrm {H}}}}^{*} {{\mathrm {H}}}{\widetilde{{{\mathrm {y}}}}}^{n+1}_+\\&= (\mathbb {1} + \gamma \Delta {{\mathrm {t}}}{{{\mathrm {H}}}}^{*} {{\mathrm {H}}}){\widetilde{{{\mathrm {y}}}}}^{n+1}_+. \end{aligned} \end{aligned}$$

We now have to work our way to (41a). First, we remark that

$$\begin{aligned} \dfrac{{\widetilde{{{\mathrm {y}}}}}^{n+1}_- - {\widetilde{{{\mathrm {y}}}}}^{n}_+}{\Delta {{\mathrm {t}}}} = \dfrac{{{\mathrm {y}}}((n+1)\Delta {{\mathrm {t}}}) - {{\mathrm {y}}}(n\Delta {{\mathrm {t}}})}{\Delta {{\mathrm {t}}}} - \dfrac{{\widehat{{{\mathrm {y}}}}}^{n+1}_- - {\widehat{{{\mathrm {y}}}}}^n_+}{\Delta {{\mathrm {t}}}}. \end{aligned}$$

(42)

Using centered Taylor expansions we can easily see that

$$\begin{aligned} \dfrac{{{\mathrm {y}}}((n+1)\Delta {{\mathrm {t}}}) - {{\mathrm {y}}}(n\Delta {{\mathrm {t}}})}{\Delta {{\mathrm {t}}}} = ({{\mathrm {A}}}+ \upeta {{\mathrm {V}}})\dfrac{{{\mathrm {y}}}((n+1)\Delta {{\mathrm {t}}}) + {{\mathrm {y}}}(n\Delta {{\mathrm {t}}})}{2} + \varepsilon ^n, \end{aligned}$$

(43)

with $\varepsilon ^n=O\bigl (\Delta {{\mathrm {t}}}^2\left\| \dddot{{{\mathrm {y}}}}\right\| _{C^3([n\Delta {{\mathrm {t}}},(n+1)\Delta {{\mathrm {t}}}],{\mathcal {Y}})}\bigr )$. Therefore, feeding equation (42) with (43) and (38a), we obtain

$$\begin{aligned} \dfrac{{\widetilde{{{\mathrm {y}}}}}^{n+1}_- - {\widetilde{{{\mathrm {y}}}}}^{n}_+}{\Delta {{\mathrm {t}}}} = ({{\mathrm {A}}}+ \upeta {{\mathrm {V}}}) \dfrac{{{\mathrm {y}}}((n+1)\Delta {{\mathrm {t}}}) - {\widehat{{{\mathrm {y}}}}}^{n+1}_-}{2} + ({{\mathrm {A}}}+ \upeta {{\mathrm {V}}})\dfrac{{{\mathrm {y}}}(n\Delta {{\mathrm {t}}}) - {\widehat{{{\mathrm {y}}}}}^{n}_+}{2}+ \varepsilon ^n, \end{aligned}$$

(44)

hence, (41a) holds. $\square $

We can now establish the energy estimate associated with (41a)–(41b)

Proposition B.2

The norm of the estimation error, namely

$$\begin{aligned} {\widetilde{{\mathcal {E}}}}^{n}_+ = \dfrac{1}{2} \left\| {\widetilde{{{\mathrm {y}}}}}^{n}_+\right\| ^2_{{\mathcal {Y}}}, \end{aligned}$$

satisfies the following estimate

$$\begin{aligned} \dfrac{{\widetilde{{\mathcal {E}}}}^{n+1}_+ - {\widetilde{{\mathcal {E}}}}^n_+}{\Delta {{\mathrm {t}}}}= & {} - \upeta \left\| \big (\sqrt{-{{\mathrm {V}}}}\big ) \dfrac{{\widetilde{{{\mathrm {y}}}}}^{n+1}_- + {\widetilde{{{\mathrm {y}}}}}^n_+}{2}\right\| ^2_{{\mathcal {Y}}} - \gamma \left\| {{\mathrm {H}}}{\widetilde{{{\mathrm {y}}}}}^{n+1}_+\right\| ^2_{{\mathcal {Z}}}\nonumber \\&- \gamma ^2\dfrac{\Delta {{\mathrm {t}}}}{2}\left\| {{{\mathrm {H}}}}^{*} {{\mathrm {H}}}{\widetilde{{{\mathrm {y}}}}}^{n+1}_+\right\| ^2_{{\mathcal {Z}}} + (\varepsilon ^n,\dfrac{{\widetilde{{{\mathrm {y}}}}}^{n+1}_- + {\widetilde{{{\mathrm {y}}}}}^n_+}{2})_{{\mathcal {Y}}}. \end{aligned}$$

(45)

Proof

Denoting ${\widetilde{{\mathcal {E}}}}^{n+1}_- = \dfrac{1}{2} \left\| {\widetilde{{{\mathrm {y}}}}}^{n+1}_-\right\| ^2_{{\mathcal {Y}}}$, we have, from system (41a)-(41b),

$$\begin{aligned} {\left\{ \begin{array}{ll} \dfrac{{\widetilde{{\mathcal {E}}}}^{n+1}_- - {\widetilde{{\mathcal {E}}}}^{n}_+}{\Delta {{\mathrm {t}}}} = - \upeta \left\| \big (\sqrt{-{{\mathrm {V}}}}\big ) \dfrac{{\widetilde{{{\mathrm {y}}}}}^{n+1}_- + {\widetilde{{{\mathrm {y}}}}}^n_+}{2}\right\| ^2_{{\mathcal {Y}}} + \Bigl (\varepsilon ^n,\dfrac{{\widetilde{{{\mathrm {y}}}}}^{n+1}_- + {\widetilde{{{\mathrm {y}}}}}^n_+}{2}\Bigr )_{{\mathcal {Y}}}\\ \dfrac{{\widetilde{{\mathcal {E}}}}^{n+1}_+ - {\widetilde{{\mathcal {E}}}}^{n+1}_-}{\Delta {{\mathrm {t}}}} = - \gamma \Bigl ({{{\mathrm {H}}}}^{*} {{\mathrm {H}}}{\widetilde{{{\mathrm {y}}}}}^{n+1}_+,\dfrac{{\widetilde{{{\mathrm {y}}}}}^{n+1}_+ + {\widetilde{{{\mathrm {y}}}}}^{n+1}_-}{2}\Bigr )_{{\mathcal {Y}}}\!\! \end{array}\right. } \end{aligned}$$

(46)

Equation (41b) leads to

$$\begin{aligned} \dfrac{{\widetilde{{\mathcal {E}}}}^{n+1}_+ - {\widetilde{{\mathcal {E}}}}^{n+1}_-}{\Delta {{\mathrm {t}}}}&= - \gamma ( {{{\mathrm {H}}}}^{*} {{\mathrm {H}}}{\widetilde{{{\mathrm {y}}}}}^{n+1}_+,{\widetilde{{{\mathrm {y}}}}}^{n+1}_+ + \frac{\gamma }{2} \Delta {{\mathrm {t}}}{{{\mathrm {H}}}}^{*} {{\mathrm {H}}}{\widetilde{{{\mathrm {y}}}}}^{n+1}_+)_{{\mathcal {Y}}} \\&= - \gamma \left\| {{\mathrm {H}}}{\widetilde{{{\mathrm {y}}}}}^{n+1}_+\right\| ^2_{{\mathcal {Z}}} - \gamma ^2\dfrac{\Delta {{\mathrm {t}}}}{2} \left\| {{{\mathrm {H}}}}^{*}{{\mathrm {H}}}{\widetilde{{{\mathrm {y}}}}}^{n+1}_+\right\| ^2_{{\mathcal {Y}}}, \end{aligned}$$

which, by regrouping both equations in (46), entails the desired estimate. $\square $

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

As for the linear case, one can derive the time-discrete dynamics satisfied by the estimation error.

Proposition B.3

Assuming that ${{\mathrm {y}}}\in C^3([0,T],{\mathcal {Y}})$, then the estimation error associated with (47a)–(47b) satisfies the discrete dynamical system

with the source term in (48a) is identical to the one in (41a), and the source term in (48b) is given by

$$\begin{aligned} \lambda ^n = O(\left\| {{\mathrm {y}}}((n+1)\Delta {{\mathrm {t}}}) - {\widehat{{{\mathrm {y}}}}}^e_+\right\| ^2_{{\mathcal {Y}}}). \end{aligned}$$

Proof

Similarly to the linear case, the prediction estimation error is

$$\begin{aligned} {\widetilde{{{\mathrm {y}}}}}^{n+1}_- = {{\mathrm {y}}}((n+1)\Delta {{\mathrm {t}}}) - {\widehat{{{\mathrm {y}}}}}^{n+1}_-. \end{aligned}$$

This entails (48b) since from equation (47b) and the linearization of ${{\mathrm {H}}}({{\mathrm {y}}}((n+1)\Delta {{\mathrm {t}}}))$ around the extrapolated trajectory ${\widehat{{{\mathrm {y}}}}}^e_+$ we have

$$\begin{aligned} \begin{aligned} {\widetilde{{{\mathrm {y}}}}}^{n+1}_-&= {{\mathrm {y}}}((n+1)\Delta {{\mathrm {t}}}) - {\widehat{{{\mathrm {y}}}}}^{n+1}_+ + \gamma \Delta {{\mathrm {t}}}{D{{\mathrm {H}}}^{e}_+}^{*}\big ({{\mathrm {z}}}^{n+1} - {{\mathrm {H}}}({\widehat{{{\mathrm {y}}}}}^{e}_+) - D{{\mathrm {H}}}^{e}_+({\widehat{{{\mathrm {y}}}}}^{n+1}_+ - {\widehat{{{\mathrm {y}}}}}^{e}_+)\big )\\&= (\mathbb {1} + \gamma \Delta {{\mathrm {t}}}{D{{\mathrm {H}}}^{e}_+}^{*} D{{\mathrm {H}}}^{e}_+ ){\widetilde{{{\mathrm {y}}}}}^{n+1}_+ + \Delta {{\mathrm {t}}}\lambda ^n. \end{aligned} \end{aligned}$$

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. $\square $

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

$$\begin{aligned} \dfrac{{\widetilde{{\mathcal {E}}}}^{n+1}_+ - {\widetilde{{\mathcal {E}}}}^n_+}{\Delta {{\mathrm {t}}}}= & {} - \upeta \left\| \big (\sqrt{-{{\mathrm {V}}}}\big ) \dfrac{{\widetilde{{{\mathrm {y}}}}}^{n+1}_- + {\widetilde{{{\mathrm {y}}}}}^n_+}{2}\right\| ^2_{{\mathcal {Y}}} \nonumber \\&- \gamma \left\| D{{\mathrm {H}}}^{e}_+ {\widetilde{{{\mathrm {y}}}}}^{n+1}_+\right\| ^2_{{\mathcal {Z}}} - \gamma ^2\dfrac{\Delta {{\mathrm {t}}}}{2}\left\| {D{{\mathrm {H}}}^{e}_+}^{*} D{{\mathrm {H}}}^{e}_+ {\widetilde{{{\mathrm {y}}}}}^{n+1}_+\right\| ^2_{{\mathcal {Z}}} \nonumber \\&+ \left( {\varepsilon ^n},{\dfrac{{\widetilde{{{\mathrm {y}}}}}^{n+1}_- + {\widetilde{{{\mathrm {y}}}}}^n_+}{2}}\right) _{{\mathcal {Y}}} + \left( { \lambda ^n},{\dfrac{{\widetilde{{{\mathrm {y}}}}}^{n+1}_+ + {\widetilde{{{\mathrm {y}}}}}^{n+1}_-}{2}}\right) _{{\mathcal {Y}}}. \end{aligned}$$

(49)

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 $\lambda ^n$ in the correction step. Hence, by following the demonstration of Proposition B.2 of Appendix B, one can obtain estimate (49). $\square $

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.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Cite this article

Imperiale, A., Chapelle, D. & Moireau, P. Sequential data assimilation for mechanical systems with complex image data: application to tagged-MRI in cardiac mechanics. Adv. Model. and Simul. in Eng. Sci. 8, 2 (2021). https://doi.org/10.1186/s40323-020-00179-w

Download citation

Received: 15 April 2020
Accepted: 13 October 2020
Published: 09 January 2021
DOI: https://doi.org/10.1186/s40323-020-00179-w

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

Abstract

Introduction

Position of the problem

Principles of sequential estimator design

Example of model formulation

Filtering data available from imaging modalities

Filtering 3D displacements from 3D grids

Filtering tagged-MR planes and grid

Extension of surface data

Tag planes

Tag grids

Filtering 2D apparent displacements

Time discretization using a prediction–correction scheme

Time discretization of the state observer

Characterization and iterative resolution of the correction step

Time discretization of the joint state and parameter observer

Results

Synthetic data generation

Generating tag planes and tag grids

Generating apparent displacements

Discrepancy measure in practice

Spectral analysis of the observer associated with (3D) tag planes

Estimation results

The 2-Regions case

The AHA-Regions case

Discussion

Towards clinical applications

Conclusions

References

Funding

Author information

Authors and Affiliations

Contributions

Corresponding author

Ethics declarations

Competing interests

Additional information

Publisher's Note

Appendices

Appendix A: The illuminating example of efficient state feedback in elastodynamics

Proposition A.1

Proof

Proposition A.2

Proof

Theorem A.3

Theorem A.4

Proof

Definition A.1

Theorem A.5

Proof

Appendix B: Analysis of the prediction–correction scheme

Linear model and linear observation operator

Proposition B.1

Proof

Proposition B.2

Proof

Linear model and nonlinear observation operator

Proposition B.3

Proof

Proposition B.4

Proof

Rights and permissions

About this article

Cite this article

Share this article

Keywords