Finite element apriori error estimate for a class of Cardiac Electric models

TL;DR

This paper derives a finite element a priori error estimate for the monodomain cardiac electric model coupled with nonlinear ionic models, demonstrating space-time convergence and providing computational results for various ionic models.

Contribution

It introduces a novel finite element error analysis for cardiac electric models with nonlinear ionic dynamics, establishing convergence rates and computational validation.

Findings

Achieves $o(h^2+k)$ space-time convergence rate.

Provides computational results for different ionic models.

Validates the theoretical error estimates through simulations.

Abstract

In this study we derive the Finite element apriori error estimate for the monodomain cardiac electric model in conjunction with the generic form for a class of nonlinear ionic models. The analysis establishes a $o(h^2+k)$ space-time convergence. Further the computational realization of the same for different reduced ionic test models is presented.

Tables

Table 1. Table 1: L 2 superscript 𝐿 2 L^{2} norm of the error and the space and time rates of convergence for different ionic models.

h	1/8	1/16	1/32	1/64	1/128
FHNM error	0.0153718	0.00418786	0.0010467	0.000261422	6.53429e-05
FHNM sroc	-	1.876	2.00037	2.0014	2.00027
FHNM troc	-	0.937999	1.00018	1.0007	1.00014
RMM error	0.0293018	0.00587156	0.00132513	0.000321176	7.93853e-05
RMM sroc	-	2.31917	2.14761	2.04469	2.01642
RMM troc	-	1.15959	1.07381	1.02235	1.00821
MSM error	0.0123688	0.0030894	0.000772382	0.000193005	4.82713e-05
MSM sroc	-	2.00131	1.99994	2.0000	2
MSM troc	-	1.00065	0.999971	1.00004	0.999998
APM error	0.0110065	0.00273321	0.000681779	0.000170299	4.25639e-05
APM sroc	-	2.00969	2.00322	2.00123	2.000037
APM troc	-	1.00484	1.00161	1.00062	1.00018

Equations

$$\frac{\partial v}{\partial t}-div(D(x)\nabla v)-I_{ion}(v,w)=I_{app}\Omega \times (0,T)$$

$$\frac{\partial w}{\partial t}-g(v,w)=0\Omega \times (0,T)$$

$$v(x,0)=v_0(x,0),w(x,0)=w_0(x,0)\Omega$$

$$n^{T}D(x)\nabla v=0\partial \Omega \times (0,T)$$

$$(v_t,{\psi }_1)-(\nabla .(D(x)\nabla v),{\psi }_1)=(I_{ion},{\psi }_1)$$

$$(w_t,{\psi }_2)=(g,{\psi }_2)$$

$$(v_t,{\psi }_1)-(D(x)\nabla v,\nabla {\psi }_1)=(I_{ion},{\psi }_1)$$

$$(w_t,{\psi }_2)=(g,{\psi }_2)$$

$$\displaystyle\Big{(}\sum_{i=1}^{N}v_{i,t}\phi_{i}(x),\phi_{i}\Big{)}+\Big{(}D(x)\nabla\sum_{i=1}^{N}v_{i}\phi_{i},\nabla\phi_{i}\Big{)}=(I_{ion},\phi_{j}),j=1,2,...,N.$$

$$M{v}_t+Av=I_{ion},M{w}_t=G,$$

$$(u_t,\psi {)}_X+(DGu,G\psi {)}_X=(F(u),\psi {)}_X$$

$$(u_{h,t},\psi {)}_X+(DG{u}_h,G\psi {)}_X=({\pi }_{h}f(u_h),\psi {)}_X,\forall \psi \in S_h\times S_h$$

$$\parallel u_h(T)-u(T){\parallel }_X\le \parallel u_{0,h}-u_0{\parallel }_X+C{h}^2\parallel u_0{\parallel }_X+C{h}^2{\int }_0^T(\parallel (F-{\pi }_{h}F)(u){\parallel }_X+\parallel u{\parallel }_X+\parallel u_t{\parallel }_X)dt$$

$$u-u_h=(u-R_{h}u)+(R_{h}u-u_h).$$

$$(D(x)G(R_{h}u-u),G\psi {)}_X=0,\forall \psi \in S_h\times S_h.$$

$$({\theta }_t,\psi {)}_X+(D(x)G\theta ,G\psi {)}_X=(F(u)-{\pi }_{h}F(u_h),\psi {)}_X-({\rho }_t,\psi {)}_X,$$

$$\frac{d}{dt}\parallel \theta {\parallel }_X\le \parallel (F-{\pi }_{h}F)(u){\parallel }_X+M\parallel u-u_h{\parallel }_X+\parallel {\rho }_t{\parallel }_X.$$

$$\parallel \theta (T){\parallel }_X\le \parallel \theta (0){\parallel }_X+{\int }_0^{T}M(\parallel \theta {\parallel }_X+\parallel \rho {\parallel }_X)+\parallel (F-{\pi }_{h}F)(u){\parallel }_X+\parallel {\rho }_t{\parallel }_{X}dt.$$

$$\parallel \rho {\parallel }_X\le C{h}^2\parallel u{\parallel }_X,\parallel {\rho }_t{\parallel }_X\le C{h}^2\parallel u_t{\parallel }_X,$$

$$\parallel \theta (0){\parallel }_X\le \parallel u_{0,h}-u_0{\parallel }_X+\parallel R_h{u}_0-u_0{\parallel }_X\le \parallel u_{0,h}-u_0{\parallel }_X+C{h}^2\parallel u_0{\parallel }_X$$

$$(\overline{{\partial }_t}{U}^n,\psi {)}_X+(D(x)G{U}^n,G\psi {)}_X=(F(U^{n-1}),\psi {)}_X,\forall \chi \text{in}{S}_h.$$

$$U^n-u^n=(U^n-{\tilde{U}}^n)+(U^n-u^n)={\theta }^n+{\rho }^n$$

$$(D(x)G({\tilde{U}}^n-u^n),G\psi {)}_X=(DG{\rho }^n,G\psi {)}_X=0$$

$$(\overline{{\partial }_t}{\theta }^n,\psi {)}_X+(D(x)G{\theta }^n,G\psi {)}_X=((\overline{{\partial }_t}{U}^n,\psi {)}_X+(D(x)G{U}^n,G\psi {)}_X-(\overline{{\partial }_t}{\tilde{U}}^n,\psi {)}_X-(D(x)G{\tilde{U}}^n,G\psi {)}_X$$

$$=(F(U^{n-1}),\psi )-(u_t^n,\psi {)}_X-(\overline{{\partial }_t}{\tilde{U}}^n-u_t^n,\psi {)}_X-(D(x)G{\tilde{U}}^n,G\psi {)}_X-(D(x)G{u}^n,G\psi {)}_X+(D(x)G{u}^n,G\psi {)}_X$$

$$=(F(U^{n-1}),\psi {)}_X-(F(u^n),\psi {)}_X-(\overline{{\partial }_t}({\tilde{U}}^n-u^n),\psi {)}_X-(D(x)G({\tilde{U}}^n-u^n),G\psi {)}_X-(\overline{{\partial }_t}{u}^n-u_t^n,\psi {)}_X$$

$$\parallel F(U^{n-1})-F(u^n){\parallel }_X\le C\parallel U^{n-1}-u^n{\parallel }_X\le C(\parallel U^{n-1}-u^{n-1}{\parallel }_X+\parallel u^{n-1}-u^n{\parallel }_X)$$

$$\le C(\parallel {\theta }^{n-1}{\parallel }_X+\parallel \rho {\parallel }_X+k\parallel \overline{{\partial }_t}{u}^n{\parallel }_X)$$

$$\frac{1}{2}\overline{{\partial }_t}{\parallel {\theta }^n\parallel }_X^2+\mu {\parallel G{\theta }^n\parallel }_X^2\le C(\parallel {\theta }^{n-1}{\parallel }_X+\parallel {\rho }^{n-1}{\parallel }_X+k\parallel \overline{{\partial }_t}{u}^n{\parallel }_X)(\parallel {\theta }^n{\parallel }_X)$$

$$+(\parallel \overline{{\partial }_t}{\rho }^n{\parallel }_X+\parallel \overline{{\partial }_t}{u}^n-u_t^n{\parallel }_X)\parallel {\theta }^n{\parallel }_X$$

$$\overline{{\partial }_t}{\parallel {\theta }^n\parallel }_X^2\le C{(\parallel {\theta }^{n-1}\parallel }_X^2+{\parallel {\rho }^{n-1}\parallel }_X^2+k{\parallel \overline{{\partial }_t}{u}^n\parallel }_X^2)+{\parallel \overline{{\partial }_t}{\rho }^n{\parallel }_X}^2+{\parallel \overline{{\partial }_t}{u}^n-u_t^n{\parallel }_X}^2+C{\parallel {\theta }^n{\parallel }_X}^2$$

Full text

Keywords: Apriori Estimate; Finite Element Method; ODE-PDE system.

1 Introduction

Cardiac electrophysiology is a branch of medicine and biology. Electrophysiological models of heart describe how the electrical potential is generated in every part of the heart. These models consist of coupled ODE-PDE system. Electrical activity at the cell level is described by these ODEs and PDEs describe the same at the tissue level. One of these models is the Bidomain model which is described by two degenerate non-linear parabolic reaction- diffusion equations which are coupled with a non-linear system of ODE ionic models. The other popular model is the, Monodomain model, which is a simplified form of the Bidomain model. It consists of a non-linear parabolic reaction-diffusion equation together with a non-linear system of ODE ionic models. While both finite difference method (FDM)and Finite element method (FEM) have been used to solve these models [9, 5, 8, 4].The space-time convergence analysis through apriori error estimation has not been reported for these models.

In this paper we derive the $L^{2}$ apriori error estimate for the FE analysis of Monodomain model with general form of ionic models. In the next section we will describe the cardiac electric models. In section (3) finite element formaulation is presented. In section (4) we derive the $L^{2}$ apriori error estimate for semi-discrete and fully discrete system. In section (5) a numerical example with space-time convergence results have been presented.

2 Cardiac Electric Model

Cardiac tissue is considered as the overlapping of the intra and extracellular continuous domains such that each point in the intracellular myocardium is also in the extracellular and the two domains are connected with continuous cellular membrane. Bidomain model(BDM) [4] is the characterization of such cardiac tissue and it consists of a set of mathematical equations which describe the propagation of intra and extracellular electric potentials in cardiac tissue.

The monodomain model(MDM) [4], simplified version of the BDM, with the capacity to provide significant information similar to that of BDM, is used to calculate the action potential. This model consists of a parabolic reaction-diffusion equation coupled with a system of ODE’s which is given as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $v$ and $w$ are the action potential and gating variables respectively. $g(v,w)$ and $I_{ion}$ is given by the models at the cell level, called ionic models.

2.1 Ionic Models

Hodgkin and Huxely [6] in 1952 gave the first mathematical model that describes accurately the action potential waveform. This model is complex in nature. There are various complex ionic models such as Luo Rudy 1 [10], Beeler Reuter ,LRd [6]. For large spatial and temporal investigation of any phenomena, various reduced ionic models FitzHugh Nagumo Model(FHNM), Roger-McCulloch Model(RMM) ,Aliev Panfilov Model(APM), Mitchell Schaeffer Model(MSM) [6, 2, 1, 3, 11] have been described to substantially provide the action potential at low cost.

3 Galerkin Finite Element Method (GFEM) formulation

Weak formulation of the system (1-4) is obtained, find $v(t)\in H^{1}(\Omega)$ [6], $w(t)\in L^{2}(\Omega)orL^{\infty}(\Omega)$[6] and $\psi_{1}\in H^{1}(\Omega)$ , $\psi_{2}\in L^{2}(\Omega)$,

[TABLE]

Integration by parts together with boundary condition gives

[TABLE]

Consider $\{T_{h}\}_{h}$ be a member of a family of quasi-uniform triangularizations of $\Omega$ with ${max}_{\tau\in T_{h}}diam\tau\leq h$ and let $S_{h}$ be the corresponding finite dimensional space of continuous functions on $\Omega$ which reduces to linear functions in each of the triangles of $T_{h}$.

Let $\{\phi_{i}\}_{i=1}^{N}$ be the basis functions of $S_{h}$ so that we can write $v(x,t)=\sum_{i=1}^{N}v_{i}(t)\phi_{i}(x)$ and $w(x,t)=\sum_{i=1}^{N}w_{i}(t)\phi_{i}(x)$ and equation (5) becomes

[TABLE]

The matrix form of the system is given as

[TABLE]

where $M_{ij}=((\phi_{i},\phi_{j}))$ , $A_{ij}=(\nabla\phi_{i},\nabla\phi_{j})$, I_{ion}=\Big{(}(I_{ion},\phi_{1}),...,(I_{ion},\phi_{N})\Big{)} and G=\Big{(}(g,\phi_{1}),...,(g,\phi_{N})\Big{)}.

4 $L^{2}$ error estimate for the semi-discrete problem

Let $\Omega$ be a plane convex domain with smooth boundary and consider the problem (1-4).

Assume that $D(x)$ is symmetric and uniformly positive definite, i.e. there exist $\alpha>0$ such that $\forall x\in R^{3},\forall\xi\in R^{3}$, $\xi^{T}D(x)\xi\geq\alpha{\mid\xi\mid}^{2}$.

Introducing the inner product $(\begin{bmatrix}\alpha_{1}&\alpha_{2}\end{bmatrix},\begin{bmatrix}\beta_{1}&\beta_{2}\end{bmatrix})_{X}:=(\alpha_{1},\beta_{1})+(\alpha_{2},\beta_{2})$, and the associated norm $\parallel.\parallel_{X}$, where $(.,.)$ is the standard $L^{2}$- inner product.

Let $u={\begin{bmatrix}v&w\end{bmatrix}}^{t}$ , $G=\begin{bmatrix}\nabla&0\\ 0&0\end{bmatrix}$ , $Gu=G\begin{bmatrix}v\\ w\end{bmatrix}=\begin{bmatrix}\nabla v\\ 0\end{bmatrix}$. So,the weak form of the problem becomes Find $u\in H^{1}(\Omega)\times L^{2}(\Omega)$ such that

[TABLE]

where $F(u)={\begin{bmatrix}I_{ion}(u)&g(u)\end{bmatrix}}^{t}$. Let $\pi_{h}\colon H^{1}(\Omega)\times L^{2}(\Omega)\longrightarrow S_{h}\times S_{h}$ be the projection and $u_{h}={\begin{bmatrix}v_{h}&w_{h}\end{bmatrix}}^{t}$ be the solution of the finite element formulation

[TABLE]

Theorem 1

Let $u$ be the solution of the problem $(\ref{mv}-\ref{in})$, and $u_{h}$ be the solution of the problem in semi-discrete case. Then, if $I_{i}on(v,w)$ and $g(v,w)$ is Lipschitz continuous in $v$ and $w$ and $D(x)$ is symmetric and positive definite, we have

[TABLE]

Proof: Decompose the error

[TABLE]

where $R_{h}u$ is the elliptic projection of the $u$ defined as,

[TABLE]

Now we will bound the $\theta=R_{h}u-u_{h}$ and $\rho=(u-R_{h}u)$ separately. In order to bound $\theta$, note that,

[TABLE]

Choose $\psi=\theta$,and applying Chauchy-schwartz and lipschitz continuity of F along with boundedness of $L^{2}$ projection, we get

[TABLE]

After integration we get,

[TABLE]

Now apply Gronwall’s lemma, and the following bounds for $\parallel\rho\parallel_{X}$ and $\parallel\rho_{t}\parallel_{X}$ taken from the elliptic theory [7]

[TABLE]

Also,

[TABLE]

Using (12-14) in 10, we arrive at the estimate 9.

4.1 $L^{2}$ error estimate for the fully discrete problem

Let $k$ be the time step, $t_{n}=nk$, and let $U^{n}$ be the approximation of $u(t_{n})$ in $S_{h}\times S_{h}$. We will use backward Euler Galerkin scheme and linearize the problem 7 by replacing $U^{n}$ by $U^{n-1}$ to obtain

[TABLE]

where $\bar{{\partial}_{t}}U^{n}=\frac{1}{k}(U^{n}-U^{n-1})$.

Theorem 2

Let $U^{n}$ and $u$ be solutions of 15 and 1, 2 respectively. Then, if $I_{ion}(v,w)$ and $g(v,w)$ are Lipschitz continuous in $v$ and $w$ and $D(x)$ is symmetric and uniformly positive definite, we have $\Arrowvert U^{n}-u(t_{n})\Arrowvert_{X}\leq C\Arrowvert u_{0,h}-u_{0}\Arrowvert_{X}+C(u)(h^{2}+k)\forall t_{n}\in\bar{J}$

Proof. $u^{n}=u(t_{n})$ ,

[TABLE]

where $\widetilde{U}^{n}$ is the elliptic projection of $u^{n}$ defined as

[TABLE]

$\rho^{n}$ will be bounded by Elliptic theory. Now we need to bound only $\theta^{n}$. For $\psi\in S_{h}\times S_{h}$,

[TABLE]

[TABLE]

[TABLE]

Now using Lipschitz continuity of $F$

[TABLE]

[TABLE]

Take $\psi=\theta^{n}$ and use the ellipticity of D and (18) we get,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Using the definition 15

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We have $T_{j}\leq C(u)(h^{2}+k)^{2}$. Hence,$\parallel\theta^{n}\parallel_{X}\leq C\parallel\theta^{0}\parallel_{X}+C(u)(h^{2}+k)$.

5 Numerical Test and Discussion

We consider the monodomain model with the following different ionic models in a square domain $[-1.25,1.25]^{2}$.

• •

FHNM: $I_{ion}=u(u-0.1)(1-u)-w,g(v,w)=u-2w$,

• •

RMM: $I_{ion}=u(u-0.1)(1-u)-vw,g(v,w)=u-2w$,

• •

APM: $I_{ion}=-ku(u-a)(u-1)-uw,$ $g(u,w)=\epsilon^{\prime}(-ku(u-1-a)-w)$, where $\epsilon^{\prime}=\epsilon_{0}+\mu_{1}w/(u+\mu_{2})$,

• •

MSM: $I_{ion}=-\frac{w}{\tau_{in}}u^{2}(u-1)-\frac{u}{\tau_{out}},$ $g(u,w)=\begin{cases}\frac{1-w}{\tau_{open}}&u\leq u_{gate},\\ \frac{-w}{\tau_{close}}&u>u_{gate}.\\ \end{cases}$

.

We solved it using linear finite elements in space and Backward Euler in time and the nonlinear terms $I_{ion}$ and $g$ are linearized by taking values at the previous time step. We compute the $L^{2}$ norm of the error and the space and time rate of convergence. Initial conditions for $v$ and $w$ are chosen to be 0.2 , 0.1 respectively for all the cases and $I_{app}=0,dt=h^{2},Dx[n]=h,DT[n]=dt.$ Space rate of convergence(sroc) = $log(L^{2}(e_{n-1})/L^{2}(e_{n}))/log(Dx[n-1])/Dx[n]$, time rate of convergence(troc) = $log(L^{2}(e_{n-1})/L^{2}(e_{n})/log(DT[n-1])/DT[n]$, where $e_{n}:n^{th}$ level error.

From Table 1 it is clear that as the grid resolution is increased not only the error in $L^{2}$ norm decreases but also the theortically predicted sroc and troc are realized in all the cases.

6 Conclusion

For the Monodomain CEM in a generic framework for ionic models an apriori error estimate under fem approach has been theortically established and computationally verified.

Acknowledgement

We would like to thank to the DST for the support through Inspire Fellowship.