A second-generation computational modeling of cardiac electrophysiology: response of action potential to ionic concentration changes and metabolic inhibition

Abstract

Background

Cardiac arrhythmias are becoming one of the major health care problem in the world, causing numerous serious disease conditions including stroke and sudden cardiac death. Furthermore, cardiac arrhythmias are intimately related to the signaling ability of cardiac cells, and are caused by signaling defects. Consequently, modeling the electrical activity of the heart, and the complex signaling models that subtend dangerous arrhythmias such as tachycardia and fibrillation, necessitates a quantitative model of action potential (AP) propagation. Yet, many electrophysiological models, which accurately reproduce dynamical characteristic of the action potential in cells, have been introduced. However, these models are very complex and are very time consuming computationally. Consequently, a large amount of research is consecrated to design models with less computational complexity.

Results

This paper is presenting a new model for analyzing the propagation of ionic concentrations and electrical potential in space and time. In this model, the transport of ions is governed by Nernst-Planck flux equation (NP), and the electrical interaction of the species is described by a new cable equation. These set of equations form a system of coupled partial nonlinear differential equations that is solved numerically. In the first we describe the mathematical model. To realize the numerical simulation of our model, we proceed by a finite element discretization and then we choose an appropriate resolution algorithm.

Conclusions

We give numerical simulations obtained for different input scenarios in the case of suicide substrate reaction which were compared to those obtained in literature. These input scenarios have been chosen so as to provide an intuitive understanding of dynamics of the model. By accessing time and space domains, it is shown that interpreting the electrical potential of cell membrane at steady state is incorrect. This model is general and applies to ions of any charge in space and time domains. The results obtained show a complete agreement with literature findings and also with the physical interpretation of the phenomenon. Furthermore, various numerical experiments are presented to confirm the accuracy, efficiency and stability of the proposed method. In particular, we show that the scheme is second-order accurate in space.

Background

Cardiac disease is the leading cause of deaths worldwide. A proportion of them is caused by rhythm irregularities of the heart, such as atrial fibrillation. In the healthy heart, the cardiac contraction is produced by softly propagating non-linear electrical waves of excitation. Any disturbance in conduction or coordination of electrical signals can result in abnormal heart rhythms, so called arrhythmias. Bradycardia, tachycardia, heart block, and atrial and ventricular fibrillation are examples of arrhythmias. The stimulation of cardiac cells is instigated by a sudden change in the electrical potential across the cell membrane due to the transmembrane flux of charged ions. The release and propagation of an electrical signal, which is ensured by controlled opening and closing of ions channels, is one of the most important functions of the cell. About fifty-two years ago, the first continuous mathematical model of cardiac cell designed to reproduce cell membrane action potentials is presented by Hodgkin and Huxley[1]. Ever since, many complex models have been developed for cardiac cells inspired by their approach. Most of these models can be classified in three sets. 1) “First generation” of ionic models which are able to reproduce basic ionic currents such as the Beeler-Ruter (BR)[2] and Luo-Rudy-I (LR-I)[3] models. 2) “Second generation” of models, which in addition to a biophysically detailed description of ion channel, pump and exchanger currents, also contain the intracellular ionic concentrations such as the DiFrancesco-Noble[4]. 3) simplified models that only contain the minimum set of phenomenological currents necessary to reproduce mesoescopic features of cell dynamics, e.g., conduction velocity (CV) restitution and action potential (AP) restitution[5, 6]. Generally, simulations based on first and second generation models are computationally demanding, however it is often desirable to designate the minimum key characteristics necessary to characterize a specific phenomenon and then proceed by using simplified models.

The purpose of this paper is to design a new computer model of cardiac action potential, which can be classified in the set of second generation models as in addition to a detailed description of ion channels, it also includes detailed description of intracellular concentrations. However, given the trend of researchers in the field, which aims to offer models with less computational complexity, we made a further simplification in the ions representation, the so called mean-field approximation of ionic solution, in which ions are not considered as microscopic discrete entities but as continuous charge densities. This leads us to a fully continuous model, Nernst-planck equations for the species concentrations and a modified cable equation for the electrical potential. This model has the advantage of besides containing a detailed biophysical description of the cardiac activity is less computationally demanding because of the simplification we made. This model is more general than those in literature of membrane transport as it extends them in three topics: 1) it’s a multidimensional model, 2) it allow accessing both time and space domains for the transport equation and also for the electrical potential equation, 3) it includes different reaction kinetics terms.

Introduction

In this paper we consider a class of models of cardiac cell signaling, which includes ions migration through biological membranes. The electrical propagation is described by a new cable equation by assuming that the membrane acts as an electrical circuit in which resistance and capacitance are arranged in a parallel circuit. The incorporated migrations exist for most living cells and some biochemical processes, where the motion of ions is supposed due to diffusion and to the effect of the electrical field. Furthermore ions undergo biochemical reactions. So the ions concentrations satisfy the Nernst-Planck equations, including a kinetic reaction terms and the potential is given by a new cable equation. The model is given by:

$$\left\{\begin{array}{l}\frac{\partial C_i}{\mathrm{\partial t}}-d_i\Delta C_i-m_i\mathit{\text{div}}(C_i\nabla \varphi )=F_i\left(C_1,\dots ,C_{N_s}\right)\text{in}{Q}_T\text{for}i=1,\dots ,\mathit{\text{Ns}}\\ \frac{\mathrm{\partial \varphi }}{\mathrm{\partial t}}-\mathrm{D\Delta \varphi }+\mathrm{D.}{F}_a\sum _{i=1}^{N_s}{z}_i{C}_i=D{\varphi }_{\mathit{\text{rest}}}-\frac{{\delta }_{\left(0,0\right)}{I}_{\mathit{\text{stim}}}}{C_m}& \text{in}{Q}_T\\ d_i\frac{\partial C_i}{\mathrm{\partial \upsilon }}+m_i{C}_i\frac{\mathrm{\partial \varphi }}{\mathrm{\partial \upsilon }}=0& \text{in}\sum _T\\ \varphi (t,x)=0& \text{in}\sum _T\\ C_i(0,x)=C_{i,0}\left(x\right)& \text{on}\Omega \\ \varphi (0,x)={\varphi }_0\left(x\right)& \text{on}\Omega \end{array}\right.$$

(1)

where Q _T =]0,T[×Ω,${\sum }_T=\left]0,T\right[\times \partial \Omega$, T>0; Ω is an open regular set of${\mathbb{R}}^2$ which represents the biological cell and ∂ Ω represents the cell membrane. Ω is supposed time independent because we are not interested in cell volume control, so we neglected variations of cell size. For each i, C _i is the concentration of the A _i species which has diffusion coefficient d _i , mobility m _i and valency z _i . ϕ is the electrical potential, D is a positive diffusion coefficient, F _a is Faraday constant, δ_(0,0) is a dirac at the point (0,0) which represents the center of the cell, C _m is the membrane capacitance, ϕ _rest is the resting potential, and F _i are reaction terms. For each i, F _i denotes the production rate of the species A _i due to all homogeneous reactions in which it is involved. We suppose that F _i depends continuously on the C _j ’s, and that d _i is a positive constant for each i.

$$\text{For all}i=1,\dots ,N_s{C}_{i,0}\in L^2\left(\Omega \right)\text{and satisfy}{C}_{i,0}\ge 0$$

(2)

In this paper we present a numerical simulation of such systems, for a large class of reaction kinetics, however, for the applications we considered a suicide substrate reaction. This article is organized in the following way. The next section is devoted to the modeling of the problem. Then, we did a finite element discretization of the mathematical model. After that we present applications, results and numerical experiments showing the accuracy, efficiency and stability of the proposed method. Finally, conclusions are drawn in the last section.

Governing equations

Modeling the electromigration of ions

Let us consider a cardiac cell which fills the bounded open set Ω of${\mathbb{R}}^N$, N≥ 1. This type of reaction within the membrane always contain electroactive ions${\left(A_i\right)}_{1\le i\le N_s}$ as one of their major components. The movement of ions is supposed to be due to diffusion and also the effect of electrical field. The mass conservation equation for the species A _i is

$$\frac{\partial C_i}{\mathrm{\partial t}}+\mathit{\text{div}}\left(J_i\right)=F_i,$$

(3)

where C _i is the concentration of species A _i , F _i denotes the production rate of A _i due to all the homogeneous reactions in which it is involved and J _i is its molar transport flux. Let’s mention that the concentration C _i can also change because of the production or destruction of the species A _i within the cell, which justify the production rate F _i that represents the net rate of increase of A _i (production-destruction). When F _i is positive, the cell is a source (leading to an increase in the total amount of A _i ), and when F _i is negative, it’s a sink. The functions F _i are called also source functions. Furthermore, migration is included along with diffusion as possible modes of transport of each A _i . The molar flux J _i then becomes

$$J_i=-d_i\nabla C_i-m_i{C}_i\nabla \mathrm{\varphi .}$$

(4)

There is a relationship between the mobility m _i and the diffusion constant d _i (see[7]), which is given by

$$m_i=\frac{d_i{z}_i{F}_a}{{\mathit{\text{RT}}}_e},$$

(5)

where z _i F _a is the charge carried by a mole of species A _i , R is the universal gas constant and T _e is the local temperature. The transport equation of each A _i becomes

$$\frac{\partial C_i}{\mathrm{\partial t}}-d_i\Delta C_i-m_i\mathit{\text{div}}(C_i\nabla \varphi )=F_i\left(C_1,\mathrm{..},C_{\mathit{\text{Ns}}}\right).$$

(6)

Modeling the electrical potential

The cell membrane can be modeled as a capacitor related in parallel with variables resistances and batteries. The capacitance is due to the phospholipid bilayer that separates the ions on the inside and the outside of the cell. The resistances and batteries represent the different ionic currents. Consequently, the electrophysiological behavior of a single cell is governed by the following differential equation[1]

$$\frac{\mathrm{\partial \varphi }}{\mathrm{\partial t}}=-\frac{I_{\mathit{\text{ion}}}+I_{\mathit{\text{stim}}}}{C_m},$$

where ϕ is the voltage, t is the time, I _ion represents the sum of the ionic currents flowing across the cell membrane, I _stim is an applied external stimulation current, and C _m is cell capacitance per unit surface area. Analogously, a 2D single cardiac cell can be modeled as a continuous system with the following partial differential equation[1]

$$\frac{\mathrm{\partial \varphi }}{\mathrm{\partial t}}=-\frac{I_{\mathit{\text{ion}}}+I_{\mathit{\text{stim}}}}{C_m}+\frac{1}{{\rho }_x{S}_x{C}_m}\frac{{\partial }^2\varphi }{\partial x^2}+\frac{1}{{\rho }_y{S}_y{C}_m}\frac{{\partial }^2\varphi }{\partial y^2},$$

where ρ _x and ρ _y are respectively the cellular resistivity in the x and y directions, S _x and S _y are respectively the surface-to-volume ratio in the x and y directions, and I _ion represents the total ionic current. In 2D, our purpose wasn’t to study the effects of anisotropy, consequently we take ρ=ρ _x =ρ _y and similarly S=S _x =S _y . By considering a “diffusion” coefficient D=1/(ρ S C _m ), it follows

$$\frac{\mathrm{\partial \varphi }}{\mathrm{\partial t}}=\mathrm{D\Delta \varphi }-\frac{I_{\mathit{\text{ion}}}+I_{\mathit{\text{stim}}}}{C_m}.$$

In the present model the ionic current depends on all ions existing in the cell, then

$$I_{\mathit{\text{ion}}}=\frac{F_a}{\mathrm{\rho S}}\sum _{i=1}^{N_s}{z}_i{C}_i-\frac{{\varphi }_{\mathit{\text{rest}}}}{\mathrm{\rho S}},$$

by specifying that the stimulation current is applied at the center of the cell (0,0), it follows that

$$\frac{\mathrm{\partial \varphi }}{\mathrm{\partial t}}-\mathrm{D\Delta \varphi }+\mathrm{D.}{F}_a\sum _{i=1}^{N_s}{z}_i{C}_i=D{\varphi }_{\mathit{\text{rest}}}-\frac{{\delta }_{(0,0)}{I}_{\mathit{\text{stim}}}\left(t\right)}{C_m}.$$

Variational formulation of the problem

In order to show the numerical formulation of the problem, let V=L²(0,T;H¹(Ω)) be the space of approximate solutions and W=H¹(Ω) be the space of tests functions. Let W^h be a finite element space of Lagrange P 1 included in W and V^h=L² (0,T;W^h) be the finite dimensional subspace of V. The Faedo-Galerkin formulation for the problem is given by, finding C_i,h∈V _h for i=1,…,N _s and ϕ _h ∈V _h such that ϕ _h =0 in ∂ Ω:

$$\left\{\begin{array}{l}\bullet \text{for every}{w}_h\in W_h\text{a.e. t}\in ]0,T[\text{and for}1\le i\le N_s,\\ \frac{d}{\mathit{\text{dt}}}{\int }_{\Omega }{C}_{i,h}{w}_h+d_i{\int }_{\Omega }\nabla C_{i,h}\nabla w_h+m_i{\int }_{\Omega }{C}_{i,h}\nabla {\varphi }_h\nabla w_h-{\int }_{\mathrm{\partial \Omega }}\left(d_i\frac{\partial C_{i,h}}{\mathrm{\partial \upsilon }}+m_i{C}_{i,h}\frac{\partial {\varphi }_h}{\mathrm{\partial \upsilon }}\right)w_h\\ ={\int }_{\Omega }{F}_i(C_{1,h},\dots ,C_{N_s,h})w_h\\ C_{i,h}(0,x)=C_{i,0}^h\left(x\right)\text{on}\Omega \\ \bullet \text{for all}{v}_h\in W_h\text{such that}{v}_h=0\text{in}\mathrm{\partial \Omega }\text{and a.e. t}\in ]0,T[\\ \frac{d}{\mathit{\text{dt}}}{\int }_{\Omega }{\varphi }_h{v}_h=-D{\int }_{\Omega }\nabla {\varphi }_h\nabla v_h-\mathrm{D.}{F}_a{\int }_{\Omega }\left(\sum _{j=1}^{N_s}{z}_j{C}_{j,h}\right)v_h+{\int }_{\Omega }\left(D{\varphi }_{\mathit{\text{rest}}}-\frac{{\delta }_{\left(0,0\right)}}{C_m}{I}_{\mathit{\text{stim}}}\left(t\right)\right)v_h\\ {\varphi }_h(0,x)={\varphi }_0^h\left(x\right)\text{on}\Omega \end{array}\right.$$

where$C_{i,0}^h,{\varphi }_0^h$ are the projections of C_i,0,ϕ₀ on W _h .

According to the boundary conditions we have

$${\int }_{\mathrm{\partial \Omega }}\left(d_i\frac{\partial C_{i,h}}{\mathrm{\partial \upsilon }}+m_i{C}_{i,h}\frac{\partial {\varphi }_h}{\mathrm{\partial \upsilon }}\right)w_h=0.$$

∀w _h ∈W _h , for 1≤i≤N _s , we define the function

$$G_i^{w_h}\left(C_{1,h},\dots ,C_{N_s,h},{\varphi }_h\right)=-d_i{\int }_{\Omega }\nabla C_{i,h}\nabla w_h-m_i{\int }_{\Omega }{C}_{i,h}\nabla {\varphi }_h\nabla w_h+{\int }_{\Omega }{F}_i\left(C_{1,h},\dots ,C_{N_s,h}\right)w_h.$$

For all v _h ∈W _h such that v _h =0 in ∂ Ω, we define the function

$$\begin{array}{l}{H}^{v_h}\left(C_{1,h},\dots ,C_{N_s,h},{\varphi }_h\right)=-D{\int }_{\Omega }\nabla {\varphi }_h\nabla v_h-\mathrm{D.}{F}_a{\int }_{\Omega }\left(\sum _{j=1}^{N_s}{z}_j{C}_{j,h}\right)v_h\\ +{\int }_{\Omega }\left(D{\varphi }_{\mathit{\text{rest}}}-\frac{{\delta }_{(0,0)}}{C_m}{I}_{\mathit{\text{stim}}}\left(t\right)\right)v_h\end{array}$$

The weak formulation of the system becomes, finding C_i,h∈V _h for i=1,…,N _s and ϕ _h ∈V _h such that ϕ _h =0 in ∂ Ω:

$$\left\{\begin{array}{l}\bullet \text{for every}{w}_h\in W_h\text{a.e. t}\in ]0,T[\text{and for}1\le i\le N_s,\\ \frac{d}{\mathit{\text{dt}}}{\int }_{\Omega }{C}_{i,h}{w}_h=G_i^{w_h}\left(C_{1,h},\dots ,C_{N_s,h},{\varphi }_h\right),\\ C_{i,h}(0,x)=C_{i,0}^h\left(x\right)\text{on}\Omega \\ \bullet \text{for all}{v}_h\in W_h\text{such that}{v}_h=0\text{in}\mathrm{\partial \Omega }\text{and a.e. t}\in ]0,T[\\ \frac{d}{\mathit{\text{dt}}}\underset{\Omega }{\int }{\varphi }_h{v}_h=H^{v_h}\left(C_{1,h},\dots ,C_{N_s,h},{\varphi }_h\right),\\ {\varphi }_h(0,x)={\varphi }_0^h\left(x\right)\text{on}\mathrm{\Omega .}\end{array}\right.$$

(7)

Channel blockers in the treatment of cardiac arrhythmias

Channel blockers (CB’s) are a type of drugs which binds to the enzyme inside the pore of a specific type of ion channel and blocks the flux of ions through it. Channel blockers are useful agents in antiarrhythmic drug therapy, especially supraventricular tachyarrhythmias[8–10]. Furthermore, there is many genetic diseases that modify and block cardiac ion channels, causing cardiac channelopathies[11]. Consequently, to model such a channel inhibition, we need to establish the transport system for suicide substrate reaction.

Suicide substrate kinetics

An enzyme system of major experimental concern; see[12, 13], is the mechanism-based inhibitor, or suicide substrate system, represented by Walsh et al.[14],

$$\begin{array}{cc}E+S\stackrel{k_1}{\underset{k_{-1}}{\rightleftharpoons }}X\stackrel{k_2}{\to }& Y\stackrel{k_3}{\to }E+P\\ Y\stackrel{k_4}{\to }{E}_i\end{array}$$

(8)

where E, S and P stand for enzyme, substrate, and product, respectively; X and Y, enzyme-substrate intermediates; E _i , inactivated enzyme; and the k’s are positive rate constants.

In this system, Y has a choice of one of two pathways, namely, to E+P with rate k₃ or to E _i with rate k₄. The ratio of these rates, k₃/k₄, is called the partition ratio and is denoted by r. Each of these pathways are supposed to be irreversible over the timescale of the reaction see[15]. S is known as a suicide substrate because it binds to the active site of an enzyme—like a substrate—but the enzyme converts it into an inhibitor which irreversibly inactivates the enzyme. Thereby, the enzyme ‘commits suicide’. In this way, a suicide substrate can specifically target an enzyme for inactivation. Furthermore, suicide substrates are particularly useful in drug administration, as they are not noxious in their common form and only the designated enzyme can convert them to their inhibitor form. For example, suicide substrates have been subject of investigation for use in the treatment of depression (monoamine oxidase inhibitors, Seiler et al.[12]), epilepsy (brain GABA transaminase inhibitors, Walsh[13]), and some tumors (ornithine decarboxylase inhibitors, Seiler et al.[12]). Suicide substrate kinetics have been studied by Waley[15] and by Tatsunami et al.[16], who had interest in the factor which determined whether the substrate was exhausted before all the enzyme was inactivated. Waley proposed it was r μ, where μ is the ratio of the initial concentration of enzyme to that of substrate, namely, e₀/s₀. Tatsunami et al., on the other hand, found the determining factor to be (1+r)μ. When (1+r)μ>1 the substrate is exhausted, while for (1+r)μ<1, all the enzyme is inactivated. When (1+r)μ=1, both occur. The interest is when e₀/s₀ is not small, which was in effect assumed since both Waley[15] and Tatsunami et al.[16] used a quasi-steady state approximation. The validity decreases for increasing values of e₀/s₀. We denote the concentrations of the reactants by

$$C_1=\left[E\right],C_2=\left[S\right],C_3=\left[X\right],C_4=\left[Y\right],C_5=\left[E_i\right],C_6=\left[P\right]$$

The law of mass action applied to (8) leads to one equation for each reactant and hence a system of nonlinear equations. We obtain the following system

$$\left\{\begin{array}{ll}\frac{\partial C_1}{\mathrm{\partial t}}-d_1\Delta C_1-m_1\mathit{\text{div}}(C_1\nabla \varphi )=-k_1{C}_1{C}_2+k_{-1}{C}_3+k_3{C}_4& \text{in}{Q}_T\\ \frac{\partial C_2}{\mathrm{\partial t}}-d_2\Delta C_2-m_2\mathit{\text{div}}(C_2\nabla \varphi )=-k_1{C}_1{C}_2+k_{-1}{C}_3& \text{in}{Q}_T\\ \frac{\partial C_3}{\mathrm{\partial t}}-d_3\Delta C_3-m_3\mathit{\text{div}}(C_3\nabla \varphi )=k_1{C}_1{C}_2-\left(k_{-1}+k_2\right)C_3& \text{in}{Q}_T\\ \frac{\partial C_4}{\mathrm{\partial t}}-d_4\Delta C_4-m_4\mathit{\text{div}}(C_4\nabla \varphi )=k_2{C}_3-\left(k_3+k_4\right)C_4& \text{in}{Q}_T\\ \frac{\partial C_5}{\mathrm{\partial t}}-d_5\Delta C_5-m_5\mathit{\text{div}}(C_5\nabla \varphi )=k_4{C}_4& \text{in}{Q}_T\\ \frac{\partial C_6}{\mathrm{\partial t}}-d_6\Delta C_6-m_6\mathit{\text{div}}(C_6\nabla \varphi )=k_3{C}_4& \text{in}{Q}_T\\ \frac{\mathrm{\partial \varphi }}{\mathrm{\partial t}}-\mathrm{D\Delta \varphi }+\mathrm{D.}{F}_a\sum _{i=1}^{N_s}{z}_i{C}_i=D{\varphi }_{\mathit{\text{rest}}}-\frac{{\delta }_{\left(0,0\right)}}{C_m}{I}_{\mathit{\text{stim}}}\left(t\right)& \text{in}{Q}_T\\ d_i\frac{\partial C_i}{\mathrm{\partial \upsilon }}+m_i{C}_i\frac{\mathrm{\partial \varphi }}{\mathrm{\partial \upsilon }}=0& \text{in}\sum _T,\text{for}i=1,2,3,\dots ,6.\\ \varphi (t,x)=0& \text{in}\sum _T\\ C_1(0,x)=e_0,C_2(0,x)=s_0,C_i(0,x)=0& \text{on}\Omega ,\text{for}i=3,\dots ,6.\\ \varphi (0,x)={\varphi }_0\left(x\right)& \text{on}\Omega \end{array}\right.$$

(9)

Numerical scheme

In this section, we present the numerical scheme for solving the problem, we used for the time marching scheme an implicit scheme for the transport equations and an explicit second order Runge-Kutta[17] scheme for the potential equation.

Time marching scheme

Our method is based on an explicit second order Runge-Kutta scheme for the potential equation and an implicit scheme for the transport equations. To this end, let us denote by$\left(C_{1,h}^{n+1},C_{2,h}^{n+1},\dots ,C_{6,h}^{n+1},{\varphi }_h^{n+1}\right)$ and$\left(C_{1,h}^n,C_{2,h}^n,\dots ,C_{6,h}^n,{\varphi }_h^n\right)$ the approximate value at time t=tⁿ⁺¹ and t=tⁿ, respectively and by δ t the time step size. Then by using (7) and the following algorithm, we determine the unknown fields.

Algorithm of resolution

We used the following algorithm to calculate ϕ _h and C_i,h.

• Initialize for i=1,…,6

  $$C_{i,h}^0=C_{i,0}^h\left(x\right),$$
  $${\varphi }_h^0={\varphi }_0^h\left(x\right)$$

• Loop over n

At step n:

• Calculate${\varphi }_h^{n+1}$ solution of

  $${\int }_{\Omega }{\varphi }_h^{n+1}{v}_h={\int }_{\Omega }\left({\varphi }_h^n+\frac{{\varphi }_h^{k_1}+{\varphi }_h^{k_2}}{2}\right)v_h.$$

  where:

  $${\int }_{\Omega }{\varphi }_h^{k_1}{v}_h=\mathrm{\delta t}·H^{v_h}\left(C_{1,h}^n,C_{2,h}^n,\dots ,C_{6,h}^n,{\varphi }_h^n\right);$$

  and

  $${\int }_{\Omega }{\varphi }_h^{k_2}{v}_h=\mathrm{\delta t}·H^{v_h}\left(C_{1,h}^n+C_{1,h}^{k_1},C_{2,h}^n+C_{2,h}^{k_1},\dots ,C_{6,h}^n+C_{6,h}^{k_1},{\varphi }_h^n+{\varphi }_h^{k_1}\right);$$

  where:

  $${\int }_{\Omega }{C}_{i,h}^{k_1}{w}_h=\mathrm{\delta t}·G_i^{w_h}\left(C_{1,h}^n,C_{2,h}^n,\dots ,C_{6,h}^n,{\varphi }_h^n\right)\forall i=1,\dots 6.$$

• Calculate$C_{1,h}^{n+1},C_{2,h}^{n+1},C_{3,h}^{n+1},C_{4,h}^{n+1},C_{5,h}^{n+1},C_{6,h}^{n+1}$ solutions of:

  initialize

  $C_{i,h}^{n+1,0}=C_{i,h}^n$ for i=1,…,6,

  Loop over k untill

  $$\sum _{i=1}^6{∥C_{i,h}^{n+1,k+1}-C_{i,h}^{n+1,k}∥}_{L_2\left(\Omega \right)}<\mathit{\text{eps}}$$
  $$\begin{array}{l}{\int }_{\Omega }\frac{C_{1,h}^{n+1,k+1}-C_{1,h}^n}{\mathrm{\delta t}}{w}_h+d_1{\int }_{\Omega }\nabla C_{1,h}^{n+1,k+1}\nabla w_h+m_1{\int }_{\Omega }{C}_{1,h}^{n+1,k+1}\nabla {\varphi }_h^{n+1}\nabla w_h\\ ={\int }_{\Omega }\left(-k_1{C}_{2,h}^{n+1,k}{C}_{1,h}^{n+1,k+1}+k_{-1}{C}_{3,h}^{n+1,k}+k_3{C}_{4,h}^{n+1,k}\right)w_h\end{array}$$
  $$\begin{array}{l}{\int }_{\Omega }\frac{C_{2,h}^{n+1,k+1}-C_{2,h}^n}{\mathrm{\delta t}}{w}_h+d_2{\int }_{\Omega }\nabla C_{2,h}^{n+1,k+1}\nabla w_h+m_2{\int }_{\Omega }{C}_{2,h}^{n+1,k+1}\nabla {\varphi }_h^{n+1}\nabla w_h\\ ={\int }_{\Omega }\left(-k_1{C}_{2,h}^{n+1,k+1}{C}_{1,h}^{n+1,k+1}+k_{-1}{C}_{3,h}^{n+1,k}\right)w_h\end{array}$$
  $$\begin{array}{l}{\int }_{\Omega }\frac{C_{3,h}^{n+1,k+1}-C_{3,h}^n}{\mathrm{\delta t}}{w}_h+d_3{\int }_{\Omega }\nabla C_{3,h}^{n+1,k+1}\nabla w_h+m_3{\int }_{\Omega }{C}_{3,h}^{n+1,k+1}\nabla {\varphi }_h^{n+1}\nabla w_h\\ ={\int }_{\Omega }\left(k_1{C}_{2,h}^{n+1,k+1}{C}_{1,h}^{n+1,k+1}-\left(k_{-1}+k_2\right)C_{3,h}^{n+1,k+1}\right)w_h\end{array}$$
  $$\begin{array}{l}{\int }_{\Omega }\frac{C_{4,h}^{n+1,k+1}-C_{4,h}^n}{\mathrm{\delta t}}{w}_h+d_4{\int }_{\Omega }\nabla C_{4,h}^{n+1,k+1}\nabla w_h+m_4{\int }_{\Omega }{C}_{4,h}^{n+1,k+1}\nabla {\varphi }_h^{n+1}\nabla w_h\\ ={\int }_{\Omega }\left(k_2{C}_{3,h}^{n+1,k+1}-\left(k_3+k_4\right)C_{4,h}^{n+1,k+1}\right)w_h\end{array}$$
  $${\int }_{\Omega }\frac{C_{5,h}^{n+1,k+1}-C_{5,h}^n}{\mathrm{\delta t}}{w}_h+d_5{\int }_{\Omega }\nabla C_{5,h}^{n+1,k+1}\nabla w_h+m_5{\int }_{\Omega }{C}_{5,h}^{n+1,k+1}\nabla {\varphi }_h^{n+1}\nabla w_h={\int }_{\Omega }{k}_4{C}_{4,h}^{n+1,k+1}{w}_h$$
  $${\int }_{\Omega }\frac{C_{6,h}^{n+1,k+1}-C_{6,h}^n}{\mathrm{\delta t}}{w}_h+d_6{\int }_{\Omega }\nabla C_{6,h}^{n+1,k+1}\nabla w_h+m_6{\int }_{\Omega }{C}_{6,h}^{n+1,k+1}\nabla {\varphi }_h^{n+1}\nabla w_h={\int }_{\Omega }{k}_3{C}_{4,h}^{n+1,k+1}{w}_h$$

Results and discussion

In this section, aiming to understand how an action potential emerges from the mathematical structure that we have developed we study the dynamics of the model for different types of input. Pulse input, time-dependent input, and sinusoidally varying amplitude input are considered in turn. These input scenarios have been chosen so as to provide an intuitive understanding of dynamics of the model. We present the behavior of the Voltage in response to a short current input, a time dependent input, and a sinusoidal current input. To describe signaling in a cell body, this one can be assumed to be an ellipse. For all the results of this section, we considered the following parameters:

For the computations cell capacitance per unit surface area is taken as C _m =2.0 μ F/c m² and surface to volume ratio is set to S=0.2 μ m^-1, following Bernus et al.[18]. Taggart et al.[19] found the velocity for conductance along the fiber direction in human myocardium 70 c m/s, which required a cellular resistivity ρ=162Ω c m for Tusscher et al.[20]. This is of the same magnitude of ρ=180Ω c m used by Bernus et al.[18] and the ρ=181Ω c m used by Jongsma and Wilders[21], and it leads to a diffusion coefficient D=1/(ρ S C _m ) of 0.00154 c m²/m s. The cell is represented by an ellipse with semi-major axis a = 2 and semi-minor axis b = 1. The diffusion coefficients of the ions are d₁=10^-3m².s^-1, d₂=2.10^-3m².s^-1, d₃=5.10^-3m².s^-1, d₄=10^-3m².s^-1, d₅=2.10^-3m².s^-1, d₆=4.10^-6m².s^-1, the reaction parameters are k₁=2 s^-1, k_-1=4 s^-1, k₂=12 s^-1, k₃=10 s^-1 and k₄=2 s^-1. The charge number of the ions are z₁=1, z₂=-1,z₃=1, z₄=1, z₅=1 and z₆=-1. The initial concentrations are e₀=0.5 μ M and s₀=0.5 μ M;ϕ _rest =0 and ϕ₀=-80 m V. The data employed for the reaction parameters and initial concentrations were taken from Burke et al.[22]. The time step of the simulation is d t=10^-5s, and T=0.006 s. The stimulus current I _stim , is the key to stimulate the system. In the heart, the excitation is ensured by the Sino-Atrial Node. Here we applied a single stimulus, which delivers a short current pulse of 1 m s and strength -200 μ A/c m², beginning at t=120×10^-5s at the center of the cell. Let’s mention that the Dirac at (0.0) can be approached by f= exp(-1000(x²+y²)).

Numerical result 1: electrophysiological validation of the model

To validate the model three criterion are considered 1) excitability 2) All or none 3) action potential morphology.

Excitability

This means that the cell is in its resting potential as long as no stimulation is applied to it. However, by using an efficient stimulus it produces AP. Figure1 shows that before time t=120.10^-5s there was no stimulation and consequently there is no action potential but at time t=120.10^-5s an AP is generated.

Figure 1

Generated AP by stimulus current at time t =120 . 10^-5s . Before time t=120.10^-5s no stimulation is applied.

All or none

This means that if the amplitude of the stimulus pulse is equivalent to the threshold, an action potential is generated, but if it’s lesser than threshold AP is not generated. Figure2 shows using a stimulus current which is less than threshold AP does not generate.

Figure 2

If the stimulus current is smaller than threshold mount, AP does not generate.

Action potential morphology

This model has characteristics of cardiac cell as it reproduces the triangular AP morphology (see Figure1) with no sustained plateaus, which is similar to the AP shape obtained with more complex models Nygren et al.[23] and cherry et al.[24]. Furthermore, the figure is similar to the action potential block by saxitoxin in the book[25].

Numerical result 2: response to a time dependent current input

In order to explore a more realistic input scenario, we stimulate the model by a time dependant input current$I_{\mathit{\text{stim}}}\left(t\right)=-200cos\left(\frac{\mathrm{\pi t}}{\mathit{\text{dt}}}\right)\mathrm{\mu A}/c{m}^2$ of 1m s duration, starting at t=120×10^-5s at the center of the cell. It can be observed that the stimulus change of sign at each time step, the results are shown in Figure3.Figure4 shows the response of electrical potential to a stimulus with a sinusoidally varying amplitude.

Figure 3

Voltage response to a time dependent current input at the center of the cell.

Figure 4

Voltage response to current input at the center of the cell with a sinusoidally varying amplitude.

Numerical result 3: visualizing current spread and its impact on the product and substrate concentrations

Figure5 plots the spread of potential through the cell. The system is stimulated by the same stimulus defined in the beginning of the section.Figure6 shows the impact of the current spread on the product concentration as it presents the spatial distribution of the product concentration from the initial time to final time.Figure7 shows the impact of the current spread on the substrate concentration as it presents the spatial distribution of the substrate concentration from initial time to final time.

Figure 5

The propagation of an electrical signal through the cell from initial time t = 0 to final time T = 0.006 s.

Figure 6

Spatial distribution of the product concentration from initial time t = 0 to final time T = 0.006 s.

Figure 7

Spatial distribution of the substrate concentration from initial time t = 0 to final time T = 0.006 s.

The convergence test

The rate of convergence of the scheme is difficult to prove analytically. However, numerical experimentation suggests that the scheme is second-order accurate in space. A quantitative estimate of the convergence error was obtained by performing a number of simulations for the same initial condition on a set of increasingly finer space meshes and time steps. The initial conditions are constants. Let T _h the mesh generation of Ω, and h(T _h )=m a x{d i a m(e _k )|e _k ∈T _h }, we take h=0.1,h=0.15 and h=0.2. For each mesh we integrate to time T with$\mathit{\text{dt}}=\frac{\mathit{\text{hT}}}{16}$. Note that as we refine the space step we also refine the time step. The error of the numerical solution was defined as$E\left(h\right)=\mathit{\text{dt}}\times \stackrel{N_s}{\sum _{i=1}}\stackrel{N_t}{\sum _{n=0}}\underset{k}{max}{∥C_{i,h}^{n,k+1}-C_{i,h}^{n,k}∥}_{L^2\left(\Omega \right)}$.

In Table1 is presented the error of convergence for different mesh sizes.

Table 1 Error of convergence for different mesh sizes

Stability and accuracy tests

Now, let us give some information about the numerical stability of our algorithm. We perform a numerical experiment with different time step dt,$\frac{\mathit{\text{dt}}}{2}$ and$\frac{\mathit{\text{dt}}}{4}$. These results suggest that the scheme is indeed stable as the solutions are quasi the same for different time steps. To illustrate, we chose to represent the electric potential, and the product concentration. In Figure8, we display snapshots of the electrical potential at time T = 0.006 s with three different time steps dt = 0.000005, dt = 0.0000025 and dt = 0.00001. We can see that the results are quasi the same at the final time T.In Figure9, we display snapshots of the product concentration at time T = 0.006 s with three different time steps dt = 0.0000025, dt = 0.000005, and dt = 0.00001. We can see that the results are quasi the same at the final time T.In Figure10, we present the evolution of the electrical potential at the center of the cell for three different time steps dt = 0.0000025, dt = 0.000005, and dt = 0.00001. We can see that the results are quasi the same.

Figure 8

Snapshots of the electrical potential at T = 0.006 s with three different time steps. The time steps are shown below each figure. (a) Time step dt = 0.00001 (b) Time step dt = 0.000005 (c) Time step dt = 0.0000025. Legend: values of concentrations through the cell.

Figure 9

Snapshots of the product concentration at T = 0.006 s with three different time steps. The time steps are shown below each figure. (a) Time step dt = 0.00001 (b) Time step dt = 0.000005 (c) Time step dt = 0.0000025. Legend: values of concentrations through the cell.

Figure 10

The evolution of the electrical potential for three different time steps. The time steps are shown below each figure. (a) Time step dt = 0.00001 (b) Time step dt = 0.000005 (c) Time step dt = 0.0000025.

Conclusion

In this paper, a new second-generation model of the set simulating cardiac action potential is proposed by using a more general mathematical model and a numerical technique based on the finite element method. The electrophysiological validation of the problem shows that the model has all characteristics of cardiac cells: excitability, all or none, and the triangular AP morphology which is similar to those obtained with more complex models Nygren et al.[23] and Cherry et al.[24]. Moreover, after comparison we can observe a complete consistency with literature findings[25]. A variety of numerical experiments were presented to confirm the accuracy, efficiency, and stability of the proposed method.