On the existence-uniqueness and computation of solution of a coupled PDE-ODE system with application to cardiac electric activity

TL;DR

This paper investigates the existence, uniqueness, and numerical computation of solutions for a coupled PDE-ODE system modeling cardiac electrical activity, combining theoretical analysis with finite element methods and simulations.

Contribution

It provides a rigorous proof of existence and uniqueness for a degenerate reaction-diffusion system in cardiac modeling, and develops a finite element based numerical scheme with simulations.

Findings

Global existence of solutions established

Uniqueness proven under specific nonlinear conditions

Numerical solutions successfully computed with FreeFem++

Abstract

In this study, we consider a system of degenerate reaction-diffusion equations, which govern the electric activity in the heart with a diffusion term modeling the potential in the surrounding tissue and the nonlinear ionic model proposed by Morris $\&$ Lecar. The global existence of a solution is established based on regularization argument using Fedo-Galerkin/Compactness approach. The uniqueness of the solution is shown based on Gronwell's Lemma upon some special treatment of nonlinear terms. The system of the continuous space-time model is first reduced to a semi-discrete time-dependent system based on finite element formulation, and then the fully discrete system is derived using the Backward Euler time stepping scheme. The numerical solution obtained using FreeFem++ are presented.

Tables

Table 1. Table 1: parameter Values used for the simulation

parameter	$v_1$	$v_2$	$v_3$	$v_4$	$v_{Ca}$	$v_K$
values	-1.2	18	-1	14.5	120	-70

Table 2. Table 2: parameter Values used for the simulation

parameter	$v_L$	$g_l$	$g_K$	$g_{Ca}$	$\lambda$	${\sigma }_l$	${\sigma }_t$
values	-50	4	8	3	1	$1.2\ast {10}^{\left(-3\right)}$	$2.5562\ast {10}^{\left(-4\right)}$

Equations

$$D_{i,e}(x)={{\sigma }_l}^{i,e}{a}_l(x){a_l}^T(x)+{{\sigma }_t}^{i,e}{a}_t(x){a_t}^T(x)+{{\sigma }_n}^{i,e}{a}_n(x){a_n}^T(x).$$

$$D_{i,e}(x)={{\sigma }_t}^{i,e}I+({{\sigma }_l}^{i,e}-{{\sigma }_t}^{i,e})a_l(x){a_l}^T(x).$$

$$C_m\frac{\partial v}{\partial t}-div(D_i(x)\nabla u_i)+I_{ion}(v,w)=I_{app}^{i}in\Omega \times (0,T)$$

$$C_m\frac{\partial v}{\partial t}-div(D_e(x)\nabla u_e)+I_{ion}(v,w)=I_{app}^{e}in\Omega \times (0,T)$$

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

$$v(x,0)=u_i(x,0)-u_e(x,0)=v_0(x,0),w(x,0)=w_0(x,0).$$

$${\int }_{\Omega }{I}_{app}^i={\int }_{\Omega }{I}_{app}^e$$

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

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

$$D(x)={\sigma }_{t}I+({\sigma }_l-{\sigma }_t)a_l(x){a_l}^T(x)\text{with}{\sigma }_{l,t}=\frac{\lambda {\sigma }_{l,t}^i}{1+\lambda }.$$

$$-I_{ion}(v,w)=\frac{1}{C_m}(g_{Ca}{m}_{\infty }(v)(v-v_{Ca})+g_{K}w(v-v_K)+g_L(v-v_L))$$

$$g(v,w)=\varphi \frac{w-w_{\infty }(v)}{{\tau }_w(v)}$$

$$C_m{\int }_{\Omega }{\partial }_{t}v{\varphi }_i+{\int }_{\Omega }{D}_i\nabla u_i\nabla {\varphi }_i+{\int }_{\Omega }{I}_{ion}(v,w){\varphi }_i={\int }_{\Omega }{I}_{app}{\varphi }_i$$

$$C_m{\int }_{\Omega }{\partial }_{t}v{\varphi }_e-{\int }_{\Omega }{D}_e\nabla u_e\nabla {\varphi }_e+{\int }_{\Omega }{I}_{ion}(v,w){\varphi }_e={\int }_{\Omega }{I}_{app}{\varphi }_e$$

$${\partial }_{t}w+g(v,w)=0$$

$$V_n=<{\{h_k\}}_{k=1}^n>$$

$$C_m{\int }_{\Omega }{\partial }_t{v}_{n}h+\frac{1}{n}{\int }_{\Omega }{\partial }_t{u}_{i,n}h+{\int }_{\Omega }{D}_i\nabla u_{i,n}\nabla h+{\int }_{\Omega }{I}_{ion}(v_n,w_n)h={\int }_{\Omega }{I}_{app}h$$

$$C_m{\int }_{\Omega }{\partial }_t{v}_{n}e-\frac{1}{n}{\int }_{\Omega }{\partial }_t{u}_{e,n}e-{\int }_{\Omega }{D}_e\nabla u_{e,n}\nabla e+{\int }_{\Omega }{I}_{ion}(v_n,w_n)e={\int }_{\Omega }{I}_{app}e$$

$${\partial }_t{w}_n+g(v_n,w_n)=0a.e.in{Q}_T$$

$$v_n(0)=v_{0,n},u_{i.n}(0)=u_{i,0,n},u_{e.n}(0)=u_{e,0,n},w_n(0)=w_{0}a.e.in\Omega$$

$$v_{0,n},u_{i,0,n},u_{e,0,n},w_{0,n}\to v_0,u_{i,0},u_{e,0},w_0\text{in}{V}^3\times L^2(\Omega ).$$

$${\int }_{\Omega }{\partial }_t{w}_n+{\int }_{\Omega }g(v_n,w_n)=0$$

$$M\left[\begin{array}{c}{c}_i^{{}^{\prime }}\\ c_e^{{}^{\prime }}\\ c_w^{{}^{\prime }}\end{array}\right]=\left[\begin{array}{c}{G}_i(t,c_i,c_e,c_w)\\ G_e(t,c_i,c_e,c_w)\\ G_w(t,c_i,c_e,c_w)\end{array}\right]$$

$$C_m{\int }_{\Omega }\frac{\partial v}{\partial t}\zeta -{\int }_{\Omega }D(x)\nabla v\nabla \zeta +{\int }_{\Omega }{I}_{ion}(v,w)\zeta =0$$

$${\int }_{\Omega }\frac{\partial w}{\partial t}\varsigma -{\int }_{\Omega }g(v,w)\varsigma =0$$

$$M\frac{\partial v_h}{\partial t}+A{v}_h+M{I}_{ion}^h(v_h,w_h)=0,\frac{\partial w_h}{\partial t}=g(v_h,w_h)$$

Full text

and

1 Introduction

An electrocardiogram examines for the problems with the electrical activity of your heart. Bidomain model[7] is used for describing the cardiac electrical activity. It consists of two PDEs coupled to a system of ODEs, describing the electrical activity of the heart. The two PDEs describes the dynamics of intracellular and extracellular potentials, whereas the ODEs, also known as an ionic model, describes the electrical behavior of the myocardium cell membrane. Involvement of different space and time scales makes it computationally expensive.

A simplified mathematical model of the cardiac tissue is the anisotropic monodomain system [7], which consist of a parabolic reaction-diffusion equation, describing the evolution of the membrane potential, coupled with an ionic model. Some of the well known ionic models are FHN, Winfree, Rogers and McCulloch, Panfilov, Phase-I Luo Rudy, Morris Lecar model, etc.

The first wellposedness of the Bidomain model with ionic model given by FHN ionic model [5] has been proved in [2]. In [3] and [9], the existence of the solution is proved for a wide class of ionic models (including Panfilov [1] and McCulloch [6]). Uniqueness, however, is achieved only for the FHN model. In [12] existence, Uniqueness and some regularity results are proved with LR1 ionic model [4].

None of those mentioned above works consider the Morris Lecar ionic model. So, This paper describes the existence of the Bidomain model with Morris Lecar as the ionic model. The main result states the existence of a weak global solution for Bidomain equation with ionic model Morris Lecar. Numerical simulation using the Finite Element Method for the monodomain model with Morris Lecar ionic model is also done here. For implemantaion FreeFem++ software is used.

In the next section, we describe the mathematical model (anisotropic Bidomain and monodomain models) and ionic model also. In section 3 we state our main existence result for Bidomain equation with Morris Lecar ionic model. In section 4 we give the proof of this result. Regularization argument and Fedo-Galerkin Technique is used to proof the result. In section 5 we provide some numerical results for the monodomain model with the Morris Lecar ionic model.

2 Mathematical model

2.1 Bidomain model

Bidomain model is a representation of the cardiac tissue as the superposition of two continuous anisotropic media, the intra (i) and extra (e) cellular media, coexisting at every point of the tissue and connected by a distributed continuous cellular membrane. This model describes the averaged intracellular and extracellular electric potentials and currents by a reaction-diffusion system of the degenerate parabolic type.

Let $\Omega\subset$ be the cardiac tissue domain. For the Bidomain characterization of the cardiac tissue, it is considered as the overlapping of the intracellular and extracellular continuous domains such that each point in the intracellular myocardium is also in the extracellular and the intracellular, and the extracellular medium, is identified by the conductivity tensors $D_{i}$ and $D_{e}$. Let $a_{l}(x),a_{t}(x),a_{n}(x)$ be orthonormal set of triplets corresponding to the structure of the cardiac tissue at a point x, where, $a_{l}$, parallel to the local fiber direction, and $a_{n}$, normal to the cardiac muscle sheet. Let ${\sigma_{l}}^{i,e}$,${\sigma_{t}}^{i,e}$, ${\sigma_{n}}^{i,e}$ are the conductivity coefficients along the corresponding directions. In general, these coefficients may depend on $x$, but in the following, we assume that they are constant, i.e., homogeneous anisotropy. Then the conductivity tensors $D_{i}$ and $D_{e}$, generally dependent on the position $x$, is given by

[TABLE]

If ${\sigma_{n}}^{i,e}$=${\sigma_{t}}^{i,e}$ we recover the axially isotropic case

[TABLE]

The bioelectric activity of cardiac cells is due to the flow $I_{ion}$ (per unit area of the membrane surface) of various ionic currents (the most important being sodium, potassium, and calcium) through the cellular membrane. Since the membrane behaves as a capacitor, the total membrane current per unit volume is given by $I_{m}=\chi(C_{m}\frac{\partial v}{\partial t}+I_{ion})$ where $v=u_{i}-u_{e}$ is the transmembrane potential, the coefficient $\chi$ is the ratio of membrane area per tissue volume, $C_{m}$ is the surface capacitance of the membrane, and $I_{ion}$ is the ionic current described later and depending on the membrane model.

Imposing the conservation of currents, we have $divJ_{i}=-divJ_{e}=I_{m}$ where $J_{i,e}=-D_{i,e}u_{i,e}$, are the intracellualr and the extracellular current densities. Therefore, in the Bidomain model, the intra and extracellular potentials $u_{i},u_{e}$ are modeled by the following reaction–diffusion system of PDEs, coupled with a system of ODEs for gating variables, descibed later. Given an applied current per unit volume $I_{app}\colon\Omega\times(0,T)\rightarrow R$ , initial conditions $v_{0}\colon\Omega\rightarrow R$ , $w_{0}\colon\Omega\rightarrow R^{M}$, find the intra and extracellular potentials $u_{i},u_{e}\colon\Omega\rightarrow R$ , the transmembrane potential $v=u_{i}-u_{e}$ and the gating variables $w_{0}\colon\Omega\rightarrow R^{M}$, such that

[TABLE]

[TABLE]

[TABLE]

We assume that the cardiac tissue is insulated, therefore homogeneous Neumann boundary conditions are assigned on $\partial\Omega\times(0,T)$ as $n^{T}D_{i,e}\nabla u_{i,e}=0$.

Initial conditions are assigned in $\Omega$ for $t=0$ as follows

[TABLE]

we then have the following compatibility condition for the system to be solvable:

[TABLE]

2.2 Simplified monodomain model

Assume that the anisotropy ratio of the two continuous media are equal, i.e. $D_{i}=\lambda D_{e}$ with $\lambda$ constant, and setting $D=\frac{\lambda D_{i}}{1+\lambda}$ and $I_{app}=\frac{\lambda I^{i}_{app}}{1+\lambda}+\frac{I^{e}_{app}}{1+\lambda}$, then the Bidomain system reduces to the anisotropic Monodomain model consisting in a parabolic reaction–diffusion equation for the transmembrane potential $v$ only which is descibed as by the following set of equations:

[TABLE]

[TABLE]

with Neumann boundary condition for $v$ and initial conditions for $v$ and $w$.

The conductivity tensor in the axial symmetric case is given by

[TABLE]

This model has been extremely used in computation because it requires substantially less computational and memory resources than the Bidomain model. Nevertheless, it is not an adequate cardiac model since it is unable to reproduce some patterns and morphology of the experimentally observed extracellular potential maps and electrograms. Therefore, unequal anisotropy ratio of the intra and extracellular media cannot be neglected.

2.3 Ionic model

Morris Lecar Model[8]

Morris Lecar model is a biological neuron model developed by Dr. Catherine Morris and Dr. Harold Lecar. A variety of oscillatory behavior of $Ca^{++}$ and $K^{+}$ conductance in the giant barnacle muscle fiber is replicated.

It consists of a two-dimensional system of non-linear ordinary differential equations. It is a simplified model version of the Hodgkin-Huxley ionic model. This system of equations qualitatively describes the complex relationship between cellular membrane potential and the ion channel activation, within the membrane: the potential depends on the activity of ion channels, the activity of ion channel depends on the voltage. Ionic current and the dynamics of gating variables are given by:

[TABLE]

[TABLE]

where

$m_{\infty}(v)=0.5[1+tanh(\frac{v-v_{1}}{v_{2}})]$

$w_{\infty}(v)=0.5[1+tanh(\frac{v-v_{3}}{v_{4}})]$

$\tau_{w}(v)=1/cosh(\frac{v-v_{3}}{2v_{4}})$

$v_{L},v_{Ca},v_{K}=$ equilibrium potential corresponding to leak , $Ca^{++},K^{+}$ conductances, respectively.

$v_{1}=$ potential at which $m_{\infty}=0.5$

$v_{2}=$ reciprocal of slope of voltage dependence of $m_{\infty}$

$v_{3}=$ potential at which $w_{\infty}=0.5$

$v_{4}=$ reciprocal of slope of voltage dependence of $w_{\infty}$

3 Existence of the bidomain model with Morris Lecar ionic model

Assumption

(A) We Assume that the conductivities of the intracellular and the extracellular spaces $D_{i}$, $D_{e}$ $\in[L_{\infty}(\Omega)]^{3\times 3}$ are symmetric and uniformly positive definite, i.e. there exist $\alpha_{i}>0,\alpha_{e}>0$ such that $\forall x\in R^{3},\forall\xi\in R^{3}$,

$\hskip 28.45274pt\xi^{T}D_{i}(x)\xi\geq\alpha_{i}{\mid\xi\mid}^{2}\hskip 14.22636pt,\hskip 14.22636pt\xi^{T}D_{e}(x)\xi\geq\alpha_{e}{\mid\xi\mid}^{2}$

Notation: $V=H^{1}(\Omega)=\{u:u\in L^{2}(\Omega),Du\in L^{2}(\Omega)\}$

Definition 1

(Weak Formulation) A weak solution of the Bidomain equation is a quadruplet of functions $(v,u_{i},u_{e},w)$ with the regularity

$v\in L^{\infty}(0,T;H^{1}(\Omega))\cap H^{1}(0,T;L^{2}(\Omega))$**

$u_{i},u_{e}\in L^{\infty}(0,T;H^{1}(\Omega)),w\in W^{1,\infty}(0,T;L^{\infty}(\Omega))$**

[TABLE]

[TABLE]

[TABLE]

for all $(\phi_{i},\phi_{e})\in H^{1}(\Omega)\times H^{1}(\Omega)$ Equations (9) and (10) holds in $D^{\prime}(0,T)$ and equation (11) holds a.e..

The next theorem provides the main result for the existence of solution for the Bidomain equation with Morris Lecar ionic model.

Theorem 1

*Let $T>0$, $I_{app}\in L^{2}(Q_{T}),D_{i},D_{e}\in[L^{\infty}(\Omega)]^{3\times 3}$ symmetric and satisfying Assumption (A), $v_{0}\in H^{1}(\Omega),w_{0}\in L^{\infty}(\Omega)$ be the given data. If $w_{0}\in L^{\infty}(\Omega)$ with a positive lower bound $r>0$, such that $r<w_{0}\leq 1$ in $\Omega$, then the problem ( 2-4) with (7), (8) and initial conditions)has a weak solution in the sense of Definition (1). *

The next section gives the proof of this theorem.

4 Proof of the theorem 3.1

The non-linear reaction-diffusion equations are degenerate in time. This issue is overcome here by adding a couple of Regularization terms, making bidomain equations parabolic. Regularization and approximation of solution are merged here. Then the resulting regularized system can be analyzed through a Fedo-Galerkin/compactness procedure and specific treatment of non-linear terms.

In section 4.1 , Regularization and Fedo-Galerkin Techniques are merged by introducing a regularized problem in finite dimension.

4.1 A regularized problem in finite dimension

Let ${\{h_{k}\}}_{k\in N^{*}}$ be a Hilbert basis of V. We assume that the basis functions are sufficiently smooth and that ${\{h_{k}\}}_{k\in N^{*}}$ is an orthonormal basis in $L^{2}(\Omega)$.

For all $n\in N^{*}$, we define the finite-dimensional space $V_{n}$ generated by ${\{h_{k}\}}^{n}_{k=1}$ i.e.

[TABLE]

Hence, we can introduce, for each $n\in N^{*}$, the following discrete problem associated with (9-11): Discrete Problem Find $(u_{i,n},u_{e,n})\in C^{1}(0,T;V_{n}\times V_{n}),w_{n}\in C^{1}(0,T;L^{\infty}(\Omega))$ such that, for $v_{n}=u_{i,n}-u_{e,n}$ and for all $(h,e)\in V_{n}\times V_{n}$ we have,

[TABLE]

[TABLE]

[TABLE]

And verifying the initial conditions

[TABLE]

The auxiliary initial conditions for $u_{i,n}$ and $u_{e,n}$ needed, are defined by introducing two arbitrary functions $u_{i,0}$ , $u_{e,0}\in V$ such that $v_{0}=u_{i,0}-u_{e,0}$ in $\Omega$. Then. for $n\in N^{*}$, we define $u_{i,0,n},u_{e,0,n},w_{0,n}$ as the orthogonal projections on $V_{n}\times V_{n}\times V_{n}$, of $u_{i,0},u_{e,0},w_{0}$.

By construction of these sequences, we have

[TABLE]

4.2 Local existence of the discretized solution

Lemma 4.1

Suppose that there exists $C$ such that

${\|u_{i,0,n}\|}_{H^{1}(\Omega)}$* + ${\|u_{e,0,n}\|}_{H^{1}(\Omega)}$ +${\|w_{0,n}\|}_{L^{2}(\Omega)}\leq C$*

For all $n\in N^{*}$ there exists a positive time $0<t_{n}<T$ which only depends on C such that the Discrete problem (12)-(15) admits a unique solution over the time interval $[0,t_{n}]$.

Proof. Since ${\{h_{k}\}}^{n}_{k=1}$ is a Hilbert Basis of $V_{n}$ and orthonormal basis of $L^{2}(\Omega)$. Since $L^{\infty}(\Omega)\subset L^{2}(\Omega)$ so $w_{0},w\in L^{2}(\Omega)$. Now we are choosing $\theta\in V_{n}$ and multiplying on both side of equation (14), we get

[TABLE]

Since ${\{h_{k}\}}^{n}_{k=1}$ is a Hilbert Basis of $V_{n}$ and orthonormal basis of $L^{2}(\Omega)$. So, now we can write,

$u_{i,n}=\sum^{n}_{l=1}c_{i,l}(t)h_{l}$ , $u_{e,n}=\sum^{n}_{l=1}c_{e,l}(t)h_{l}$ , $w_{n}=\sum^{n}_{l=1}c_{w,l}(t)h_{l}$ ,

$u_{i,0,n}=\sum^{n}_{l=1}c^{0}_{i,l}h_{l}$ , $u_{e,0,n}=\sum^{n}_{l=1}c^{0}_{e,l}h_{l}$ , $w_{0,n}=\sum^{n}_{l=1}c^{0}_{w,l}h_{l}$

Notation

$c_{i}=\{c_{i,l}\}^{n}_{l=1}$ , $c_{e}=\{c_{e,l}\}^{n}_{l=1}$ , $c_{w}=\{c_{w,l}\}^{n}_{l=1}$

$c^{0}_{i}=\{c^{0}_{i,l}\}^{n}_{l=1}$ , $c^{0}_{e}=\{c^{0}_{e,l}\}^{n}_{l=1}$ , $c^{0}_{w}=\{c^{0}_{w,l}\}^{n}_{l=1}$

Then the system of equations (12),(13),(16) is equivalent to the following non-linear system of ODEs

[TABLE]

$\begin{bmatrix}c_{i}(0)\\ c_{e}(0)\\ c_{w}(0)\end{bmatrix}$ = $\begin{bmatrix}c^{0}_{i}\\ c^{0}_{e}\\ c^{0}_{w}\end{bmatrix}$

Here the matrix $M\in R^{3n\times 3n}$ is given by

$M$= $\begin{bmatrix}(C_{m}+\frac{1}{n})M_{V}&\vdots-C_{m}M_{V}&\vdots&0\\ \cdots\hskip 14.22636pt\vdots&\cdots&\cdots\\ -C_{m}M_{V}&(C_{m}+\frac{1}{n})M_{V}&\vdots&0\\ \cdots&\cdots&\cdots\\ 0\hskip 28.45274pt\vdots&0\hskip 28.45274pt\vdots&M_{V}\end{bmatrix}$

with $M_{V}\in R^{n\times n}$

and $M_{V}=(\int_{\Omega}h_{k}h_{l})_{1\leq k,l\leq n}$

and right-hand side of (17) is given by

$G_{i}(t,c_{i},c_{e},c_{w})=-\int_{\Omega}D_{i}\nabla u_{i,n}\nabla h_{k}-\int_{\Omega}I_{ion}(v_{n},w_{n})h_{k}+\int_{\Omega}I_{app}h_{k}$

for all $1\leq k\leq n$,

$G_{i}(t,c_{i},c_{e},c_{w})=-\int_{\Omega}D_{i}\nabla u_{i,n}\nabla h_{k}-\int_{\Omega}I_{ion}(v_{n},w_{n})h_{k}+\int_{\Omega}I_{app}h_{k}$

for all $1\leq k\leq n$,

$G_{w}(t,c_{i},c_{e},c_{w})=-\int_{\Omega}g(v_{n},w_{n})h_{k}$

for all $1\leq k\leq n$,

Lemma 4.2

For all $n\in N^{*}$, the matrix M is positive definite.

According to Lemma (4.2) the mass matrix M is positive definite and hence invertible and, on the other hand, the RHS of (17) is a $C^{1}$ function with respect to the arguments $c_{i},c_{e},c_{w}$. So, By using the Cauchy-Lipschitz theorem, existence of the local solution of the ODE system (17) follows.

4.3 Energy estimates

Lemma 4.3

Let $u_{i,0},u_{e,0}\in H^{1}(\Omega)$ and $w_{0}\in L^{\infty}(\Omega)$ with $r<w_{0}\leq 1$ , there exists a positive constant $w_{min}$ independent of $T^{\prime}$ such that a solution $(u_{i,n},u_{e,n},w_{n})$ of the discrete problem defined on $[0,T^{\prime}]$ for $T^{\prime}>0$ satisfies

$\|v_{n}\|_{L^{\infty}(0,t;L^{2}(\Omega))}+\frac{1}{\sqrt{n}}(\|u_{i,n}\|_{L^{\infty}(0,t;L^{2}(\Omega))}+\|u_{e,n}\|_{L^{\infty}(0,t;L^{2}(\Omega))})+\|\nabla u_{i,n}\|_{L^{2}(Q_{t})}+\|\nabla u_{e,n}\|_{L^{2}(Q_{t})}\leq c$**

$\|\partial_{t}v_{n}\|_{L^{2}(Q_{t})}+\|v_{n}\|_{L^{\infty}(0,t;H^{1}(\Omega))}+\frac{1}{\sqrt{n}}(\|\partial_{t}u_{i,n}\|_{L^{2}(Q_{t})}+\|\partial_{t}u_{e,n}\|_{L^{2}(Q_{t})})+\|\nabla u_{i,n}\|_{L^{\infty}(0,t;L^{2}(\Omega))}+\|\nabla u_{e,n}\|_{L^{\infty}(0,t;L^{2}(\Omega))}\leq c$**

and, for all $t\in[0,T^{\prime}]$

$\|w_{n}\|_{W^{1,\infty}(0,t;L^{\infty}(\Omega))}\leq c$* , $w_{min}\leq w_{n}\leq 1$ in $Q_{T}$.*

**Proof.**From equation (14) it follows that

$\partial_{t}w_{n}=-g(v_{n},w_{n})$

$\partial_{t}w_{n}=-\phi\frac{w_{n}-w_{\infty}(v_{n})}{\tau(v_{n})}$ where $0\leq w_{\infty}\leq 1$

$\partial_{t}w_{n}\geq-\phi\frac{w_{n}}{\tau(v_{n})}$

$w_{n}\geq w_{0}exp(-\frac{T\phi}{\tau(v_{n})})$

again

$\partial_{t}w_{n}\leq-\phi\frac{w_{n}-1}{\tau(v_{n})}=\phi\frac{1-w_{n}}{\tau(v_{n})}$

$w_{n}\leq 1-(1-w_{0})exp(-\frac{T\phi}{\tau(v_{n})})$

using $r<w_{0}\leq 1$, we obtain that

$w_{min}=rexp(-\frac{T\phi}{\tau(v_{n})})\leq w_{n}\leq 1$, a.e. in $Q_{T}$.

On the other hand,

$-\frac{\phi}{\tau_{w}}\leq\partial_{t}w_{n}\leq\frac{\phi}{\tau_{w}}$

$\|w_{n}\|_{W^{1,\infty}(0,t;L^{\infty}(\Omega))}\leq\frac{\phi}{\tau_{w}}$

For first energy estimate take $h=u_{i,n}$ and $e=u_{e,n}$ in equation (12) and (13) and then subtract, which yields

$\frac{1}{2}\frac{d}{dt}[C_{m}\|v_{n}\|^{2}_{L^{2}(\Omega)}+\frac{1}{n}(\|u_{i,n}\|^{2}_{L^{2}(\Omega)}+\|u_{e,n}\|^{2}_{L^{2}(\Omega)})]+\alpha_{i}\|\nabla u_{i,n}\|^{2}_{L^{2}(\Omega)}$

$+\alpha_{e}\|\nabla u_{e,n}\|^{2}_{L^{2}(\Omega)}+\int_{\Omega}I_{ion}(v_{n},w_{n})v_{n}\leq\int_{\Omega}I_{app}v_{n}$

here,

$-I_{ion}(v,w)=\frac{1}{C_{m}}(g_{Ca}m_{\infty}(v)(v-v_{Ca})+g_{K}w(v-v_{K})+g_{L}(v-v_{L}))$

$I_{ion}v=-\alpha v^{2}-\beta v^{2}w-\gamma v-\delta vw$

where, $\alpha=\frac{1}{C_{m}}(g_{Ca}m_{\infty}(v)+g_{L})$ , $\beta=\frac{1}{C_{m}}g_{K}$ ,

$\gamma=\frac{1}{C_{m}}(v_{Ca}g_{Ca}m_{\infty}(v)+g_{L}v_{L})$ , $\delta=-\frac{1}{C_{m}}g_{K}v_{K}$

So, $I_{ion}v\geq-\alpha^{{}^{\prime}}\mid v\mid^{2}-\gamma^{{}^{\prime}}$ (because $\alpha$ and $\gamma$ are bounded)

where $\alpha^{{}^{\prime}}$ and $\gamma^{{}^{\prime}}$ are positive constants.

Therefore,

$\frac{1}{2}\frac{d}{dt}[C_{m}\|v_{n}\|^{2}_{L^{2}(\Omega)}+\frac{1}{n}(\|u_{i,n}\|^{2}_{L^{2}(\Omega)}+\|u_{e,n}\|^{2}_{L^{2}(\Omega)})]+\alpha_{i}\|\nabla u_{i,n}\|^{2}_{L^{2}(\Omega)}$

$+\alpha_{e}\|\nabla u_{e,n}\|^{2}_{L^{2}(\Omega)}\leq(\alpha^{{}^{\prime}}+\frac{1}{2})\|v_{n}\|^{2}_{L^{2}(\Omega)}+\frac{1}{2})\|I_{app}\|^{2}_{L^{2}(\Omega)}+\gamma^{{}^{\prime}}\mid\Omega\mid$

Therefore, integrating over $(0,t)$ with $t\in[0,T^{\prime}]$, and after that applying Gronwall Lemma and using the fact that,

$\frac{1}{n}({\|u_{i,0,n}\|}_{L^{2}(\Omega)}+{\|u_{e,0,n}\|}_{L^{2}(\Omega)})$ +${\|v_{0,n}\|}^{2}_{L^{2}(\Omega)}$ is uniformly bounded with respect to $n$, we will get first estimate.

For the estimate of the time derivatives, we take $h=\partial_{t}u_{i,n}$ and $e=\partial_{t}u_{e,n}$ in equation (12) and (13) and subtract then integrating over $(0,t)$ with $t\in[0,T^{\prime}]$, we obtain

$\frac{1}{4}C_{m}\|\partial_{t}v_{n}\|^{2}_{L^{2}(Q_{t})}+\frac{1}{n}(\|\partial_{t}u_{i,n}\|^{2}_{L^{2}(Q_{t})}+\|\partial_{t}u_{e,n}\|^{2}_{L^{2}(Q_{t})})+\frac{\alpha_{i}}{2}\|\nabla u_{i,n}\|_{L^{2}(\Omega)}+\|\nabla u_{e,n}\|_{L^{2}(\Omega)}$

$\leq c(\|\nabla u_{i,0,n}\|_{L^{2}()\Omega})+\|\nabla u_{e,0,n}\|_{L^{2}()\Omega})+\frac{1}{2}\|I_{app}\|_{L^{2}(Q_{t})}+\gamma^{{}^{\prime}}T+\alpha^{{}^{\prime}}\|v_{n}\|^{2}_{L^{2}(Q_{t})}+\int^{t}_{0}\int_{\Omega}(\beta v_{n}+\delta)w_{n}\partial_{t}v_{n}$

The last term is solved as $\int^{t}_{0}\int_{\Omega}(\beta v_{n}+\delta)w_{n}\partial_{t}v_{n}\leq\int^{t}_{0}\int_{\Omega}(\beta v_{n}+\delta)\partial_{t}v_{n}\leq\|(\beta v_{n}+\delta)\|_{L^{2}(Q_{t})}\|\partial_{t}v_{n}\|_{L^{2}(Q_{t})}$.

So, by inserting this inequality in the last expression and using the previous estimate we obtain the second estimate.

4.4 Global existence of solution

Energy estimates allows us to extend the existence time of our discrete solution $u_{i,n},u_{e,n},w_{n}$. According to Lemma (4.3),the solution satisfies, for all $t\in[0,T^{\prime}]$, where $T^{\prime}$ is the existence time

${\|u_{i,n}(t)\|}_{H^{1}(\Omega)}$ + ${\|u_{e,n}(t)\|}_{H^{1}(\Omega)}$ +${\|w_{n}(t)\|}_{L^{2}(\Omega)}\leq c$

After applying Lemma (4.1) iteratively, we obtain the existence of solution upto an arbitrary time $T$.

We now want to pass to the limit when n goes to infinity. Let us multiply equation (12) (13) by a function $\eta\in D(0,T)$ and integrate between $0\&T$. For all $k\leq n$, we have

$C_{m}\int^{T}_{0}\int_{\Omega}\eta\partial_{t}v_{n}h_{k}+\frac{1}{n}\int^{T}_{0}\int_{\Omega}\eta\partial_{t}u_{i,n}h_{k}+\int^{T}_{0}\int_{\Omega}\eta D_{i}\nabla u_{i,n}\nabla h_{k}$

$+\int^{T}_{0}\int_{\Omega}\eta I_{ion}(v_{n},w_{n})h_{k}=\int^{T}_{0}\int_{\Omega}\eta I_{app}h$

$C_{m}\int^{T}_{0}\int_{\Omega}\eta\partial_{t}v_{n}h_{k}-\frac{1}{n}\int^{T}_{0}\int_{\Omega}\eta\partial_{t}u_{e,n}h_{k}-\int^{T}_{0}\int_{\Omega}\eta D_{e}\nabla u_{e,n}\nabla h_{k}$

$+\int^{T}_{0}\int_{\Omega}\eta I_{ion}(v_{n},w_{n})h_{k}=\int^{T}_{0}\int_{\Omega}\eta I_{app}h$

From Lemma (4.3) it follows that there exists $u_{i},u_{e}\in L^{\infty}(0,T;H^{1}(\Omega)),v\in L^{\infty}(0,T;H^{1}(\Omega))\cap H^{1}(0,T;L^{2}(\Omega)),w\in L^{\infty}(Q_{T})$ such that

$u_{i,e,n}\rightarrow u_{i,e}$ in $L^{\infty}(0,T;H^{1}(\Omega))$ weak $*$

$v_{n}\rightarrow v$ in $L^{\infty}(0,T;H^{1}(\Omega))$ weak $*$

$v_{n}\rightarrow v$ in $H^{1}(0,T;L^{2}(\Omega))$ weak

$w_{n}\rightarrow$ w in $L^{\infty}(Q_{T})$ weak $*$

According to Lemma (4.3), we also conclude that $\frac{1}{\sqrt{n}}u_{i,n}$ and $\frac{1}{\sqrt{n}}u_{e,n}$ are bounded in $L^{\infty}(0,T;L^{2}(\Omega))$. So, for all $k\in N^{*}$ and $\alpha\in D(0,T)$, we have

$lim_{n\to\infty}\frac{1}{n}\int^{T}_{0}\int_{\Omega}\eta\partial_{t}u_{i,n}h_{k}=0$ , $lim_{n\to\infty}\frac{1}{n}\int^{T}_{0}\int_{\Omega}\eta\partial_{t}u_{e,n}h_{k}=0$

Since $v_{n}$ is bounded in $L^{2}(0,T;H^{1}(\Omega))\cap H^{1}(0,T;L^{2}(\Omega))$, so $v_{n}$ will be bounded in $H^{1}(Q_{T})$. Hence by compact embedding of $H^{1}(Q_{T})$ in $L^{3}(Q_{T})$, the sequence $\{v_{n}\}$ strongly converges to $v$ in $L^{3}(Q_{T})$,

$lim_{n\to\infty}\int^{T}_{0}\int_{\Omega}\eta(\alpha(v_{n})v_{n}+\gamma(v_{n}))h_{k}=\int^{T}_{0}\int_{\Omega}\eta(\alpha(v)v+\gamma(v))h_{k}$

Since sequence $\{w_{n}\}$ is bounded in $L^{\infty}(Q_{T})$ and $v_{n}$ strongly converges to $v$ in $L^{2}(Q_{T})$, we have

$lim_{n\to\infty}\int^{T}_{0}\int_{\Omega}\eta(\beta v_{n}+\delta)w_{n}h_{k}=\int^{T}_{0}\int_{\Omega}\eta(\beta v+\delta)wh_{k}$

Thus, $lim_{n\to\infty}\int^{T}_{0}\int_{\Omega}\eta I_{iom}(v_{n},w_{n})h_{k}=\int^{T}_{0}\int_{\Omega}\eta I_{ion}(v,w)h_{k}$

Since

$w_{n}\rightarrow$ w in $L^{\infty}(Q_{T})$ weak $*$

4.5 Uniqueness of the weak solution

Lemma 4.4

Assume that the first partial derivatives of $I_{ion}(v,w)$ and $g(v,w)$ are bounded and that $(v_{1},u_{i,1},u_{e,1},w_{1}),(v_{2},u_{i,2},u_{e,2},w_{2})$ are two weak solutions of our problem corresponding, respectively , to the initial data $(v_{1,0},w_{1,0})$ and $(v_{2,0},w_{2,0})$ and right-hand sides $I_{app,1}$ and $I_{app,1}$. For all $t\in(0,T)$, there holds

$\|v_{1}(t)-v_{2}(t)\|^{2}_{L^{2}(\Omega)}+\|w_{1}(t)-w_{2}(t)\|^{2}_{L^{2}(\Omega)}$**

$\leq exp(K_{1}t)K_{2}(\|v_{1,0}-v_{2,0}\|^{2}_{L^{2}(\Omega)}+\|w_{1,0}-w_{2,0}\|^{2}_{L^{2}(\Omega)}+\|I_{app,1}-I_{app,2}\|^{2}_{L^{2}(Q_{t})})$**

Proof: Ref[9].

Remark: This result also provides a stability estimate with respect to the initial condition.

5 Numerical method

We consider Monodomain model with Morris Lecar ionic model and consider $I_{app}=0$. Our monodomain system is equivalent to finding $v\in H^{1}(\Omega)$ and $w\in L^{\infty}(\Omega)$ such that

[TABLE]

[TABLE]

for all $\zeta\in H^{1}(\Omega),\varsigma\in L^{2}(\Omega)$. The system is discretized in space using finite element method and in time using Backward Euler method.

5.1 Finite element discretization in space

We consider the square domain $\Omega=[a,b]^{2}$. The domain$\Omega$ is discretized by introducing a structured quasi-uniform grid of triangular $P1$ elements(which is denoted by $\tau_{h}$). So, FEM approximation for the domain $\Omega$ is, $\Omega=\cup_{E\in\tau_{h}}E$.

The associated finite element space

$V_{h}=\{\zeta_{h}\in V:\zeta_{h}\text{is continuous in}\Omega:{\zeta_{h}}_{\mid E}\in P_{1}({E}),\forall E\in\tau_{h}\}$

A semi-discrete form is obtained by applying a standard Galerkin procedure and choosing a finite element basis $\{\zeta_{i}\}\in V_{h}$. Let $I^{h}_{ion}$ be the finite element approximation of $I_{ion}$.

In the monodomain model, the finite element approximation $v_{h}$ of transmembrane potential $v$ and $w_{h}$ of gating variable $w$ are the solution of

[TABLE]

where

$M=(m_{rs}),m_{rs}=\sum_{E}\int_{E}\zeta_{r}\zeta_{s}dx$

$A=(a_{rs}),a_{rs}=\sum_{E}\int_{E}{\nabla\zeta_{r}}^{T}D(x)\nabla\zeta_{s}dx$

Numerical quadrature in 2-dimension is used in order to compute these integrals. Now this ODE system is discretized by Backward Euler and the implementation is done using FreeFem++ library functions[10].

5.2 Numerical results and discussion

To numerically simulate the transmembrane potential ($v$) in a cardiac tissue as modeled by the coupled PDE-ODE system (5-6) representing the monodomain model[7] for cardiac tissue together with the ionic model proposed by Morris and Lecar[8] the required parameter is chosen as given in[7] and[8] and these details are provided in the table 1 and table 2. In all the numerical simulations cardiac tissue has been represented by a square domain $\Omega=[-1.25,1.25]^{2}$. To carry out the numerical simulations one has to choose an appropriate grid system so that numerical solution is computed to the acceptable accuracy. Such a goal has been achieved here through grid validation test. Three different grid systems consisting of (a) 121, (b) 441 and (c) 1681 degrees of freedom (dofs) ( or (a) 200, (b) 800 and (c) 3200 number of linear triangular elements respectively) have considered. Transmembrane potential ($v$) obtained using these three grid systems are compared at different points of the domain. In all these cases only a marginal variation (less than $0.5\%$) in “$v$” is noticed as one moves from the grid system (a) to the grid system (c). As a sample in Fig. 1 the transmembrane potential corresponding to (x,y)=(0,0) with time are compared. Clearly, the grid system with 1681 dofs is more than adequate for the current set of simulations. Hence all the simulations are carried out using 3200 linear triangular elements with 1681 dofs. In Fig. 2 the temporal variation of transmembrane potential corresponding to the following five different representative points, chosen from the four different quadrants, of the domain are presented. The initial transmembrane potential at these points follow the IC setup, and they gradually evolve to the same steady state with $v\cong 2$, indicating that the intracellular and extracellular potentials have reached a state which no more favor calcium and Potassium ion migration and thereby maintain stable ionic concentration in the absence of any depolarization phase. It is noticed that in about 300-time steps the transmembrane potential corresponding to all these points already reach the steady-state and on further time marching only a marginal variation in v is noticed. So the entire domain is nearly re-polarized in about 300-time steps where each time step corresponds to 0.1-millisecond size. It is also to be noted lesser the initial transmembrane potential than its steady state value, more rapid is the growth in $v$, due to larger driving force (potential difference), and hence all points irrespective of its starting value would reach a steady state in about 300 time steps. This also indicates that smaller the $v$ than 2, stronger may be the migration of calcium and potassium ions. Now to trace the transmembrane potential in the entire domain in the form of isochrones or iso-potential plots for the entire domain corresponding to six different time instances covering the initial state to the steady state situation are presented in Fig. 3. From the legend values for $v$ it is clear that the min-max differences in $v$ reach to a nearly zero state in about 400 time-steps each of 0.1 milliseconds.

6 Conclusions

We prove the Existence uniqueness for the solution of the bidomain model with Morris Lecar ionic model using Galerkin and Compactness approach. Also, Finite Element Computations based on Monodomain model depict the success and effectiveness of Morris Lecar ionic model in enabling the visualization the Cardiac Electric Activity in cardiac tissue.

7 Acknowledgement

We would like to thank the DST for the support through Inspire Fellowship. Also thankful to FreeFem++ open source developers[11] for making it available for the researchers.