State Estimation for Electrophysiology Models with Front Level-Set Data Using Topological Gradient Derivations

, and the primary unknowns are the extracellular potential and the transmembrane potential . Hence, the intracellular potential reads . For all with $\int _{\mathcal S} u_e \, dS= 0$ such that, for all $(\phi ,\psi )$ with $\int _{\mathcal S} \psi \, dS= 0$ ,

$\begin{aligned} {\left\{ \begin{array}{ll} A_m \displaystyle \int _{\mathcal S} \Bigl ( C_m \frac{\partial u}{\partial t} + I_{\text {ion}}(u) \Bigr ) \phi \, dS + \displaystyle \int _{\mathcal S} \Bigl ({\varvec{\sigma }}_{i} \cdot \bigl (\varvec{\nabla }u+ \varvec{\nabla }u_e\bigr )\Bigr )\cdot \varvec{\nabla }\phi \, dS \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad = A_m \displaystyle \int _{\mathcal S} I_{\text {app}} \phi \, dS, \\ \displaystyle \int _{\mathcal S} \Bigl (({\varvec{\sigma }}_{i} + {\varvec{\sigma }}_{e}) \cdot \varvec{\nabla }u_e \Bigr ) \cdot \varvec{\nabla } \psi \, dS + \displaystyle \int _{\mathcal S} \Bigl ({\varvec{\sigma }}_{i} \cdot \varvec{\nabla }u\Bigr ) \cdot \varvec{\nabla } \psi \, dS = 0. \\ \end{array}\right. } \end{aligned}$

(1)

Here, the positive constant

denotes the ratio of membrane area per unit volume,

is the membrane capacitance per unit surface, $I_{\text {ion}}(u)$ a reaction term representing the ionic current across the membrane and also depending on local ionic variables satisfying additional ODEs – in our case we use the ionic model of [8] – and $I_{\text {app}}$ a given prescribed stimulus current, when applicable. We define the intra- and extra-cellular diffusion tensors ${\varvec{\sigma }}_{i}$ and ${\varvec{\sigma }}_{e}$ by

$\begin{aligned} {\varvec{\sigma }}_{i, e} = \sigma _{i,e}^t {\varvec{I}} + (\sigma _{i,e}^l - \sigma _{i,e}^t) \bigl [ I_0(\theta ) \varvec{\tau }_0 \otimes \varvec{\tau }_0 + J_0(\theta ) \varvec{\tau }_0^\perp \otimes \varvec{\tau }_0^\perp \bigr ], \end{aligned}$

(2)

where ${\varvec{I}}$ denotes the identity tensor in the tangential plane – also sometimes called the surface metric tensor – $\varvec{\tau }_0$ is a unit vector associated with the local fiber direction on the atria midsurface, and $\varvec{\tau }_0^\perp$ such that $(\varvec{\tau }_0,\varvec{\tau }_0^\perp )$ gives an orthonormal basis of the tangential plane. The functions $I_0(\theta ) = \frac{1}{2} + \frac{1}{4 \theta } \sin (2 \theta )$ and $J_0(\theta ) = 1 - I_0(\theta )$ represent the effect of an angular variation $2\theta$ of the fiber direction across the wall. A typical physiological simulation of the model is presented in Fig. 1 in a healthy case, with the parameters given in Table 1. For details on the modeling formulation and parameter calibration we refer to [4, 7], and also to [16] where this model was used in the atria for numerical simulations of complete realistic electrocardiograms.

Table 1.

Conductivity parameters (all in $\text {S.cm}^{-1}$ ) and maximal conductance $g_{Na}$ in the different atrial areas (all in $\text {nS}.\text {pF}^{-1}$ ) with RT = regular tissue, PM = pectinate muscles, CT = crista terminalis, BB = Bachman’s bundle, FO = fossa ovalis

$\sigma _e^{t}$	$\sigma _e^{l}$	$\sigma _i^{t}$	$\sigma _i^{l}$	$g_{Na}$ (RT)	$g_{Na}$ (PM)	$g_{Na}$ (CT)	$g_{Na}$ (BB)	$g_{Na}$ (FO)
$9.0\, 10^{-4}$	$2.5\, 10^{-3}$	$2.5\, 10^{-4}$	$2.5\, 10^{-3}$	7.8	11.7	31.2	46.8	3.9

Fig. 1.

Atrial electrical depolarization and corresponding synthetic front data

2.2 Data of Interest

We assume in this work that the patient-specific depolarization front is measured, as is the case when isochrones are available. From a mathematical standpoint, the measurement procedure can be modeled by considering that, for a particular solution of (1) denoted by $(\breve{u},\breve{u}_e)$ and associated with patient-specific parameters and initial conditions, we have at our disposal the time evolution of the front

$\begin{aligned} \varGamma _{\breve{u}}(t) = \{\varvec{x} \in \mathcal S, \, \breve{u}(\varvec{x},t) = c_{\text {th}}\}, \end{aligned}$

(3)

defining $c_{\text {th}}$ as a threshold value characterizing the front, and the already traveled-through region is given by

$$$\begin{aligned} \mathcal S^{\text {in}}_{\breve{u}}(t) = \{\varvec{x} \in \mathcal S, \, \breve{u}(\varvec{x},t) > c_{\text {th}}\}, \end{aligned}$$” src=”/wp-content/uploads/2016/09/A339585_1_En_46_Chapter_Equ4.gif”></DIV></DIV><br /> <DIV class=EquationNumber>(4)</DIV></DIV>up to some measurement noise. These data can be represented as a sequence of images – denoted by <SPAN id=IEq41 class=InlineEquation><IMG alt=$$z(t)$$ src=$ – that essentially take two different values inside and outside the traveled-through region. Our objective is to use the image sequence

to reconstruct the target solution $(\breve{u},\breve{u}_e)$ , in a context where the initial conditions and some physical parameters are uncertain.

3 Estimation Methodology

3.1 Sequential Estimation Principles

We consider a general dynamical system $\dot{y} = A(y,\theta ,t)$ , where y is the state variable, A the model operator, and $\theta$ some parameters of interest. In this abstract setting, we consider a specific target trajectory $$$\{ \breve{y}(t), t>0 \}$$” src=”/wp-content/uploads/2016/09/A339585_1_En_46_Chapter_IEq46.gif”></SPAN> solution of the model with initial condition <SPAN id=IEq47 class=InlineEquation><IMG alt=$ , where $y_\diamond$ is a known a priori whereas $\breve{\zeta }_y$ is unknown, and assuming the same type of decomposition for the parameters $\theta = \theta _\diamond + \breve{\zeta }_\theta$ . We further assume that we have at our disposal some indirect measurements of the target trajectory represented by the observation variable z(t). Our estimation problem consists in reconstructing the solution $$$\{ \breve{y}(t), t>0 \}$$” src=”/wp-content/uploads/2016/09/A339585_1_En_46_Chapter_IEq51.gif”></SPAN> and possibly identifying the parameters <SPAN id=IEq52 class=InlineEquation><IMG alt=$ from the data z(t).

As a prerequisite to any estimation strategy, we must be able to define – at each time – a similarity/discrepancy measure $\mathcal {D}(y,z)$ between the data z and the state variable y. When $\mathcal {D}(y,z)$ vanishes, the state is exactly compatible with the data. By contrast, when $\mathcal {D}(y,z)$ is non-zero, the data indicate that $y(t) \ne \breve{y}(t)$ .

To achieve our estimation objective, we adopt a so-called sequential strategy where we define an observer system – also known as sequential estimator system – as a new dynamical system of the form

$\begin{aligned} {\left\{ \begin{array}{ll} \dot{\hat{y}} = A(\hat{y},\hat{\theta }) + B_y(\mathcal {D},\nabla _y \mathcal {D}), \quad y(0) = y_\diamond , \\ \dot{\hat{\theta }} = B_\theta (\mathcal {D},\nabla _y \mathcal {D}), \quad \hat{\theta }(0) = \theta _\diamond . \end{array}\right. } \end{aligned}$