 Research
 Open Access
 Published:
A secondgeneration computational modeling of cardiac electrophysiology: response of action potential to ionic concentration changes and metabolic inhibition
Theoretical Biology and Medical Modellingvolume 11, Article number: 46 (2014)
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 NernstPlanck 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 secondorder 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 nonlinear 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 fiftytwo 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 BeelerRuter (BR)[2] and LuoRudyI (LRI)[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 DiFrancescoNoble[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 meanfield 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, Nernstplanck 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 NernstPlanck equations, including a kinetic reaction terms and the potential is given by a new cable equation. The model is given by:
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.
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
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 } (productiondestruction). 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
There is a relationship between the mobility m_{ i } and the diffusion constant d_{ i } (see[7]), which is given by
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
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]
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]
where ρ_{ x } and ρ_{ y } are respectively the cellular resistivity in the x and y directions, S_{ x } and S_{ y } are respectively the surfacetovolume 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
In the present model the ionic current depends on all ions existing in the cell, then
by specifying that the stimulation current is applied at the center of the cell (0,0), it follows that
Variational formulation of the problem
In order to show the numerical formulation of the problem, let V=L^{2}(0,T;H^{1}(Ω)) be the space of approximate solutions and W=H^{1}(Ω) 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^{2} (0,T;W^{h}) be the finite dimensional subspace of V. The FaedoGalerkin 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 ∂ Ω:
where${C}_{i,0}^{h},{\varphi}_{0}^{h}$ are the projections of C_{i,0},ϕ_{0} on W_{ h }.
According to the boundary conditions we have
∀w_{ h }∈W_{ h }, for 1≤i≤N_{ s }, we define the function
For all v_{ h }∈W_{ h } such that v_{ h }=0 in ∂ Ω, we define the function
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 ∂ Ω:
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 mechanismbased inhibitor, or suicide substrate system, represented by Walsh et al.[14],
where E, S and P stand for enzyme, substrate, and product, respectively; X and Y, enzymesubstrate 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_{3} or to E_{ i } with rate k_{4}. The ratio of these rates, k_{3}/k_{4}, 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_{0}/s_{0}. 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_{0}/s_{0} is not small, which was in effect assumed since both Waley[15] and Tatsunami et al.[16] used a quasisteady state approximation. The validity decreases for increasing values of e_{0}/s_{0}. We denote the concentrations of the reactants by
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
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 RungeKutta[17] scheme for the potential equation.
Time marching scheme
Our method is based on an explicit second order RungeKutta 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^{n+1} and t=t^{n}, 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)\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{2em}{0ex}}$$ 
Loop over n
At step n:

Calculate${\varphi}_{h}^{n+1}$ solution of
$${\int}_{\Omega}{\varphi}_{h}^{n+1}{v}_{h}={\int}_{\Omega}\phantom{\rule{0.3em}{0ex}}\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}\xb7{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}\xb7{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}\xb7{G}_{i}^{{w}_{h}}\left({C}_{1,h}^{n},{C}_{2,h}^{n},\dots ,{C}_{6,h}^{n},{\varphi}_{h}^{n}\right)\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\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}{\u2225{C}_{i,h}^{n+1,k+1}{C}_{i,h}^{n+1,k}\u2225}_{{L}_{2}\left(\Omega \right)}<\mathit{\text{eps}}$$$$\begin{array}{l}\phantom{\rule{15.0pt}{0ex}}{\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}\\ \phantom{\rule{8em}{0ex}}\phantom{\rule{5em}{0ex}}={\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}$$$$\phantom{\rule{15.0pt}{0ex}}\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}\\ \phantom{\rule{8em}{0ex}}\phantom{\rule{5em}{0ex}}={\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}$$$$\phantom{\rule{15.0pt}{0ex}}\begin{array}{l}\phantom{\rule{1em}{0ex}}{\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}\\ \phantom{\rule{8em}{0ex}}\phantom{\rule{5em}{0ex}}={\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}$$$$\phantom{\rule{15.0pt}{0ex}}\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}\\ \phantom{\rule{8em}{0ex}}\phantom{\rule{5em}{0ex}}={\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}\phantom{\rule{1em}{0ex}}\right){w}_{h}\phantom{\rule{1em}{0ex}}\end{array}$$$$\phantom{\rule{15.0pt}{0ex}}{\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}$$$$\phantom{\rule{15.0pt}{0ex}}{\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, timedependent 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^{2} 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^{2}/m s. The cell is represented by an ellipse with semimajor axis a = 2 and semiminor axis b = 1. The diffusion coefficients of the ions are d_{1}=10^{3}m^{2}.s^{1}, d_{2}=2.10^{3}m^{2}.s^{1}, d_{3}=5.10^{3}m^{2}.s^{1}, d_{4}=10^{3}m^{2}.s^{1}, d_{5}=2.10^{3}m^{2}.s^{1}, d_{6}=4.10^{6}m^{2}.s^{1}, the reaction parameters are k_{1}=2 s^{1}, k_{1}=4 s^{1}, k_{2}=12 s^{1}, k_{3}=10 s^{1} and k_{4}=2 s^{1}. The charge number of the ions are z_{1}=1, z_{2}=1,z_{3}=1, z_{4}=1, z_{5}=1 and z_{6}=1. The initial concentrations are e_{0}=0.5 μ M and s_{0}=0.5 μ M;ϕ_{ rest }=0 and ϕ_{0}=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^{5}s, 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 SinoAtrial Node. Here we applied a single stimulus, which delivers a short current pulse of 1 m s and strength 200 μ A/c m^{2}, beginning at t=120×10^{5}s at the center of the cell. Let’s mention that the Dirac at (0.0) can be approached by f= exp(1000(x^{2}+y^{2})).
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^{5}s there was no stimulation and consequently there is no action potential but at time t=120.10^{5}s an AP is generated.
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.
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\phantom{\rule{0.3em}{0ex}}\left(\frac{\mathrm{\pi t}}{\mathit{\text{dt}}}\right)\mathrm{\mu A}/c{m}^{2}$ of 1m s duration, starting at t=120×10^{5}s 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.
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.
The convergence test
The rate of convergence of the scheme is difficult to prove analytically. However, numerical experimentation suggests that the scheme is secondorder 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}{\u2225{C}_{i,h}^{n,k+1}{C}_{i,h}^{n,k}\u2225}_{{L}^{2}\left(\Omega \right)}$.
In Table1 is presented the 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.
Conclusion
In this paper, a new secondgeneration 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.
References
 1.
Hodgkin A, Huxley A:A quantitative description of membrane current and its application to conduction and excitation nerve. J Physiol. 1952, 117: 500540.
 2.
Beeler G, Reuter H:Reconstruction of the ap of ventricular myocardial fiber. J Physiol. 1977, 268: 177219.
 3.
Luo C, Rudy Y:A model of the ventricular cardiac action potential. Circ Res. 1991, 68: 15011526. 10.1161/01.RES.68.6.1501.
 4.
Difrancesco D, Noble D:A model of cardiac electrical activity incorporating ionic pumps and concentraion changes. Phil Trans Roy Soc London B. 1985, 307: 353398. 10.1098/rstb.1985.0001.
 5.
Fenton F:Theoretical investigation of spiral and scroll wave instabilities underlying cardiac fibrillation. PhD thesis, Northeastern Univ. Boston University; 1999,
 6.
Fenton F, Karma A:Vortex dynamics in threedimensional continuous myocardium with fiber rotation: filament instablility and fibrillation. Chaos. 1998, 8: 2047. 10.1063/1.166311.
 7.
Keener J, Sneyd J: Mathematical Physiology. 1998, New York: Springer
 8.
Haines DE, DiMarco JP:Current therapy for supraventricular tachycardia. Curr Probl Cardiol. 1992, 17: 41110.1016/01462806(92)90011C.
 9.
Singh BN, Ellrodt G, Peter CT:Verapamil: a review of its pharmacological properties and therapeutic use. Drugs. 1978, 15: 16910.2165/0000349519781503000001.
 10.
Singh BN, Hecht HS, Nademanee K, Chew CY:Electrophysiologic and hemodynamic effects of slowchannel blocking drugs. Prog Cardiovasc Dis. 1982, 25: 10310.1016/00330620(82)900238.
 11.
Keating MT, Sanguinetti MC:Pathophysiology of ion channel mutations. Curr Opin Genet Dev. 1996, 6: 326333. 10.1016/S0959437X(96)800104.
 12.
Seiler N, Jung MJ, Kochweser J: EnzymeActivated Irreversible Inhibitors. 1978, Oxford: Elsevier Holland
 13.
Walsh CT:Suicide substrates, mechanismbased enzyme inactivators: recent developments. Annu Rev Biochem. 1984, 53: 493535. 10.1146/annurev.bi.53.070184.002425.
 14.
Walsh CT, Cromartie T, Marcotte P, Spencer R:Suicide substrates for flavprotein enzymes. Methods Enzymol. 1978, 53: 437448.
 15.
Waley SG:Kinetics of suicide substrates. Biocheml J. 1980, 53: 771773.
 16.
Tatsunami S, Yago N, Hosoe M:Kinetics of suicide substrates. steadystate treatments and computeraided exact solutions. Biochem Biophys Acta. 1981, 662: 226235.
 17.
Demailly JP: Analyse Numérique et Equations Différentielles. 1991, Grenoble: Press Universitaires
 18.
Bernus O, Wilders R, Zemlin CW, Verschelde H, Paniflov AV:A computationally efficient electrophysiological model of human ventricules. Am J Physiol Heart Circ Physiol. 2002, 282: 22962308.
 19.
Taggart P, Sutton PMI, Opthof T, Coronel R, Trimlett R, Pugsley W, Kallis P:Inhomogeneous transmural conduction during early ischemia in patients with coronary artery disease. J Mol Cell Cardiol. 2000, 32: 621639. 10.1006/jmcc.2000.1105.
 20.
ten Tusscher KH, Noble D, Noble PJ, Panfilov AV:A model for human ventricular tissue. Am J Physiol Heart Circ Physiol. 2004, 286: 15731589.
 21.
Jongsma HJ, Wilders R:Gap junctions in cardiovascular disease. Circ Res. 2000, 86: 11931197. 10.1161/01.RES.86.12.1193.
 22.
Burke MA, Maini PK, Murray JD:On the kinetics of suicide substrates. Biophysl Chem. 1990, Jeffries Wyman Anniversary Volume (37): 8190.
 23.
Nygren A, Fiset C, Firek L, Clark JW, Lindblad DS, Clark RB, Giles WR:Mathematical model of an adult human atrial cell. The role of k+ currents in repolarization. Circ Res. 1998, 82: 6381. 10.1161/01.RES.82.1.63.
 24.
Cherry EM, Hastings HM, Evans SJ:Dynamics of human atrial cell models: Restitution, memory, and intracellular calcium dynamics in single cells. Prog Biophys Mol Biol. 2008, 98: 2437. 10.1016/j.pbiomolbio.2008.05.002.
 25.
DuBois LM: Action Potential: Biophysical and Cellular Context, Initiation, Phases and Propagation. 2010, New York: Nova Science Publishers, Inc.
Acknowledgements
We are grateful to the anonymous referee for the corrections and useful suggestions that have improved this article.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
All authors contributed to writing and improving the paper and read and approved the final manuscript.
Nour Eddine Alaa, Hamid Lefraich contributed equally to this work.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Cardiac action potential
 Cable equation
 Reactiondiffusion system
 Electromigration
 Nonlinear coupled system
 Finite element method
 NernstPlanck equation
 Numerical analysis
 Substrate suicide
 Computational simulation