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

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.

http://www.tbiomed.com/content/11/1/46 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

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 http://www.tbiomed.com/content/11/1/46 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: is an open regular set of 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.
For all i = 1, . . . , N s C i,0 ∈ L 2 ( ) and satisfy C i,0 ≥ 0 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.

Modeling the electromigration of ions
Let us consider a cardiac cell which fills the bounded open set of R N , N ≥ 1. This type of reaction within the membrane always contain electroactive ions (A i ) 1≤i≤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 http://www.tbiomed.com/content/11/1/46 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 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] ∂φ ∂t 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] ∂φ ∂t 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" 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 P1 included in W and V h = L 2 0, T; W h be the finite dimensional subspace of V . The Faedo-Galerkin formulation for the problem is given by, finding The weak formulation of the system becomes, finding

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][9][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], 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 http://www.tbiomed.com/content/11/1/46 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 quasi-steady 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 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 C n+1 1,h , C n+1 2,h , . . . , C n+1 6,h , φ n+1 h and C n 1,h , C n 2,h , . . . , C n 6,h , φ n h 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 .

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 http://www.tbiomed.com/content/11/1/46 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/cm 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 cm/s, which required a cellular resistivity ρ = 162 cm for Tusscher et al. [20]. This is of the same magnitude of ρ = 180 cm used by Bernus et al. [18] and the ρ = 181 cm used by Jongsma and Wilders [21], and it leads to a diffu-  [22]. The time step of the simulation is dt = 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 Sino-Atrial Node. Here we applied a single stimulus, which delivers a short current pulse of 1 ms and strength −200 μA/cm 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. Figure 1 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. Figure 2 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 Figure 1) 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 stim (t) = −200 cos πt dt μA/cm 2 of 1ms 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 Figure 3.   Figure 4 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 Figure 5 plots the spread of potential through the cell. The system is stimulated by the same stimulus defined in the beginning of the section. Figure 6 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. Figure 7 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 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 ) = max{diam(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 dt = 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 In Table 1 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, dt 2 and dt 4 . These results sug- gest 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 Figure 8, 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 Figure 9, 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 Figure 10, 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 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.