 Research
 Open Access
 Published:
Cardiac dynamics: a simplified model for action potential propagation
Theoretical Biology and Medical Modellingvolume 9, Article number: 50 (2012)
Abstract
This paper analyzes a new semiphysiological ionic model, used recently to study reexitations and reentry in cardiac tissue [I.R. Cantalapiedra et al, PRE 82 011907 (2010)]. The aim of the model is to reproduce action potencial morphologies and restitution curves obtained, either from experimental data, or from more complex electrophysiological models. The model divides all ion currents into four groups according to their function, thus resulting into fastslow and inwardoutward currents. We show that this simplified model is flexible enough as to accurately capture the electrical properties of cardiac myocytes, having the advantage of being less computational demanding than detailed electrophysiological models. Under some conditions, it has been shown to be amenable to mathematical analysis. The model reproduces the action potential (AP) change with stimulation rate observed both experimentally and in realistic models of healthy human and guinea pig myocytes (TNNP and LRd models, respectively). When simulated in a cable it also gives the right dependence of the conduction velocity (CV) with stimulation rate. Besides reproducing correctly these restitution properties, it also gives a good fit for the morphology of the AP, including the notch typical of phase 1. Finally, we perform simulations in a realistic geometric model of the rabbit’s ventricles, finding a good qualitative agreement in AP propagation and the ECG. Thus, this simplified model represents an alternative to more complex models when studying instabilities in wave propagation.
Introduction
Cardiac action potentials (APs) are electrical signals that trigger the synchronous contraction of the heart. As such, the regular propagation of the action potential is necessary for ensuring a correct heart functioning. APs are produced as a result of ion currents that cross the cell membrane, producing a net depolarization or repolarization of the membrane as different currents are invoked in response to the transmembrane voltage changes. The currents are produced by the movement of individual ions through ion channels, which are specialized poreforming proteins that span the cell membrane, forming a pathway for ions to cross it. Each type of channel is highly selective for a specific type of ion. The most common intracellular ion concentrations considered in cardiac models are calcium, sodium, and potassium. Channelopaties due to mutations that modify the ion channel function can perturb the form of the action potential, sometimes leading to cardiac dysfunctions or altered AP propagation.
Following the pioneering mathematical description of the action potential by Hodgkin and Huxley[1] in neuronal cells, the first cardiac models considered, in much the same way, the changes in the transmembrane potential produced by a sum of ionic gated currents. As more refined experimental data of the different currents and their dynamics became available, more complex models of the AP for specific cardiac cells have been proposed[2]. The quest of new models has been especially noticeable for several animal species, as rabbit[3], guinea pig[4], rat[5], or dog[6]; and, most recently, for humans[7–10]. Nowadays, one can find very detailed models in the literature where the description of ionic channel gating is given in terms of Markov processes[11], which permit to link specific mutations to the model parameters. A known drawback of the Markov formulation is the increasing complexity of the models and consequently the very high computational cost associated to the simulation of few beats in a realistic heart geometry. Clearly, a highly realistic ionic model is ideal as a test bench, allowing to safely test multiple hypotheses and making predictions without incurring any cardiac risk for the patient, or to check conditions not easily reproducible in experiments. But, despite their more realistic description, these models are often difficult to analyze, let alone amenable for deep mathematical analysis.
In order to tackle the problem of modeling AP propagation in the computer, simplified models have also been proposed, from the very simple but unphysiological[12–14], to the socalled semiphysiological models[15–19], somewhere in between the former and the extremely detailed physiological models mentioned earlier. These semiphysiological models preserve particular sets of properties of the action potential, as for instance, the dependence of the AP duration (APD, time during which the voltage is above a certain threshold, characterizing the duration of the excited state) or the conduction velocity (CV, propagation speed of the AP pulse) on the stimulation period, known as restitution curves. In this respect, these can be viewed as mesoscopic models, that bridge the gap between dynamics at the molecular level (ion channel gating) and whole heart description, being still manageable for both computation and theoretical analysis.
A widely used semiphysiological model is the one by Fenton and Karma[15], a threevariable model of the cardiac action potential. This model uses three transmembrane currents, i.e., fast inward, slow inward, and slow outward, which resemble the set of physiological sodium, calcium, and potassium currents, respectively. In the FentonKarma model, the variables are the transmembrane potential and two gating variables that are used to regulate the inactivation of the fast inward and slow inward currents, respectively. Despite its simplicity, the model can reproduce fairly the action potential duration (APD) and conduction velocity (CV) restitution curves of more complex models or experiments. However, it fails to reproduce accurately the AP morphology. This is a strong drawback of the FentonKarma model. In fact, the model was originally created in order to describe the propagation and properties of scroll waves. The underlying idea is that wave behavior is determined by the restitution properties of the system. Thus, it has also been useful in the study of cardiac alternans (a beat to beat modification in the duration of the AP, at fast pacing rates), where the onset and subsequent evolution of alternans has been related to the shape of the restitution curves. However, more recent results have shown that models with the same restitution properties, but different AP morphologies may give different onsets for alternans[17]. The reason is that the dynamics of alternans is also affected by electrotonic currents, that depend crucially on the form of the action potential[20]. Furthermore, electrotonic currents are also a determinant factor for the occurrence of reexcitation and phase 2 reentry in cases where the duration of the action potential is nonhomogeneous in the heart tissue[21]. Typical examples are Brugada[22], or long QT[23] syndromes, that may be caused by fast or delayed repolarization, respectively.
With the aim of providing models that fit better the AP morphology, a modification of the threevariable FentonKarma model has been proposed in[19]. This new model includes an additional gate that modulates the slow inward current, in a fashion that resembles the effect of the fast outward potassium current (I_{ to }) in physiological models. This allows to reproduce the notch in phase 1 of the AP typical of epicardial cells. In contrast to[19], in the present paper we specifically include the effect of the I_{ to }. Following the same idea as in the FentonKarma model, we divide the currents not by their carrier ion, but by their function: inward or outward currents, and among these, slow and fast currents, to a total of four. With this assumption, we show that it is possible to reproduce all the important characteristics of AP propagation and morphology. Here, we specifically compare the results with more realistic models of electrical activity in human myocytes[8] and guinea pig heart ventricles[4], and with available experimental data[24], including action potential shape for different pacing rates, and different types of myocytes, i.e. epicardial; endocardial; and midmyocardial cells.
The present model has been recently used to study reexcitations in tissue presenting large dispersion of repolarization[25]. There, it was shown that reexcitation is due to a slow pulse, induced by the Ltype calcium current, that propagates from the region of long APDs to the region of short APDs until it encounters newly excitable tissue. An interesting advantage of studying this mechanism using this simplified model was the possibility to obtain the characteristics of this slow pulse analytically.
The paper is organized as follows: in the next section we introduce the equations of the model. Then, the AP morphology is studied, obtaining the proper parameters matching with experiments and two detailed models: LuoRudy (LRd) model for guineapig[4], and ten Tusscher et al. (TNNP) model of human ventricular cells[8]. A discussion of the different currents and gates follows, as well as a more detailed comparison with other semiphysiological models. In the next section, we perform simulations of AP propagation in a threedimensional model of the ventricles. In the concluding section, we provide some future lines of research opened by this paper.
Mathematical formulation of the simplified model
The cellular transmembrane potential$\stackrel{~}{V}$ (physical units are mV) satisfies the following equation:
where ∇ is the standard nabla operator,$D=\sigma /\left(\stackrel{~}{S}{C}_{m}\right)$ is the anisotropic diffusion tensor (units: cm^{2} ms^{−1}), σ is the electrical conductivity (units: mS cm^{−1}),$\stackrel{~}{S}\equiv S/\text{Vol}$ is the cell surface to volume ratio (units: cm^{−1}) and C_{ m } is the membrane capacitance per unit area (≈1 μ F cm^{−2}). The simulations in a onedimensional cable are carried out by replacing the spatial derivative term$\nabla \xb7(D\nabla \stackrel{~}{V})$ by its onedimensional counterpart$D{\partial}_{\mathit{\text{xx}}}^{2}\stackrel{~}{V}$ in Equation (1). In the threedimensional simulations, one must take into account the inhomogeneous and anisotropic behavior of the diffusion parameter[26]. All computer simulations have been performed using a simple Euler forward algorithm (for time integration) and finite spatial discretization in order to ensure current conservation up to double precision arithmetic. The time step is fixed to dt = 0. 01 ms. The uniform spatial discretization mesh used in the simulations was set to dx = 0. 02 cm for 1D and 2D, and dx = 0. 025 cm for 3D simulations. The current I_{ ion }is the sum of the ionic currents flowing across the cell membrane including the ones exchanged by the pumps through active transport (units are μ A cm^{−2}) and I_{ stim } is an applied external stimulation current. To stimulate the cell we apply it during 1 ms with an amplitude of 1.5 times the threshold stimulation value.
Realistic models include ion currents corresponding to all the existing ion channels experimentally found in the cell membrane. In the present simplified model, we decompose the total membrane current into only four components, i.e., I_{ ion }= I_{ fi } + I_{ si } + I_{ so } + I_{ to }, where the sum contains a fast inward current I_{ fi }(Na^{+} current plus the fast part of the Ca^{2 +} current); a slow inward Ca^{2+} current, I_{ si }; a slow outward timeindependent K^{+} current, I_{ so }; and a fast transient outward K^{+} current, I_{ to }. The model does not consider the intra and extracellular concentrations of different ions, and therefore the effect of the pumps is implicitly included in the previous currents. At this point, let us remind that in the FentonKarma model, the three variables of the model are the membrane potential V(x t), the inactivation gate h(t) of the fast inward current, and the inactivation gate f(t) of the slow inward current; steadystate activation is assumed for both of these currents. In the present paper, we supplement the FentonKarma model variables with two additional dynamical variables associated with the activation and inactivation of the I_{ to }current. The latter current being the essential novelty with respect to the three variable model of Fenton and Karma[15]. Let us insist that it is crucial to take into consideration this fast transient outward current to obtain accurately the characteristic notch in the phase 1 of the AP, most noticeable in epicardial myocytes. This current was experimentally studied in human hearts by Nabauer et al.[27] and Li et al.[24] and has been suggested to contribute significantly to regional electrophysiological heterogeneity in myocardial cells and tissue of several animal species. In particular, it is well known that the human ventricle shows substantial transmural heterogeneity in AP morphology related to transient outward properties. Furthermore, this current is sometimes related to the occurrence of phase2 reentry[28, 29].
In the rest of the paper, in order to alleviate the mathematical notation, it is convenient to rewrite Eq. (1) in dimensionless form by defining the membrane dimensionless potential$V=(\stackrel{~}{V}{\stackrel{~}{V}}_{\mathit{\text{res}}})/\mathrm{\Delta}\stackrel{~}{V}$, which varies roughly between zero and one and where${\stackrel{~}{V}}_{\mathit{\text{res}}}\simeq 85$ mV is the resting potential of the polarized membrane and$\mathrm{\Delta}\stackrel{~}{V}=100$ mV. One also defines the scaled currents${J}_{\mathit{\text{fi}}}={I}_{\mathit{\text{fi}}}/\left[{C}_{m}\mathrm{\Delta}\stackrel{~}{V}\right]$ (and similar expressions for the others currents: J_{ so }; J_{ to } ; J_{ si }, which have all dimension of inverse time and units are ms^{−1}). The explicit expressions for these currents are given by:
where V_{ fi }and V_{ to }are parameters indicating the reversal potentials of sodium and potassium respectively. The four gate variables obey the classical saturation kinetics:
where upper dots denote differentiation with respect to time. The expressions for the steadystate values of the gates appearing in the current formulations, and the time constants of exponential functions with which they converge, are:
where Θ(.) is the Heaviside function, V_{ r }= 0. 6, V_{ to }= 0 and the rest of the parameters of our model are given in Table1 for human experimental epi, endo, midmyocardium, and TNNP and LRd numerical models. Note that in the above expressions a + subscript refers to the time associated for a given gate to converge to one and, conversely, a − subscript refers to the time associated for a given gate to converge to zero. The determination of the numerical values of the parameters will be explained in the next section.
Results
Comparison with experimental results and detailed electrophysiological models
One of the main goals of the present model is to reproduce the AP morphology measured in experiments, or obtained with more complex mathematical models, for different types of myocytes and different pacing rates. Hence, the present model is an extension of the FentonKarma model, but with the ability of reproducing the AP morphology and, particularly, the notch in phase 1 faithfully.
We have validated our model in two ways: firstly, we have fitted experimental APs for human endo, mid, and epicardial ventricular cells choosing the appropriate model parameters. Secondly, we have also fitted the AP morphologies of two of the most currently used theoretical models in the literature: ten Tusscher et al.[8] for human ventricular myocytes (TNNP), and LuoRudy[4, 30] for guinea pig (LRd), but also widely used in other contexts. We have fitted the parameters of our model such that the AP morphologies, as well as the correct conduction velocities, at different pacing rates were obtained. The model parameters are determined by using a standard nonlinear optimization routine from Matlab based on the Levenberg Marquardt algorithm, that minimizes the square of the difference between the model action potential and the objective template. One subtlety of the optimization process used here is that we have imposed a slightly larger weight in the zone corresponding to the notch of the AP that is mostly influenced by the J_{ to }current. The codes used for the fit, as well as more information about the model can be found in[31].
Comparison with experimental human cell ventricular morphologies
As a first use of our model, we show in Figure1 the comparison between experimental data of AP morphologies and the present model. We have fitted the model to the experiments of Li et al.[24], that were done at pacing rates of 1 and 2 Hz. As it is clear from Figure1, the model reproduces well the deep notch observed in epicardial cells; the more accentuated action potential duration of Mcell; as well as the variation of the action potential duration with the stimulation rate in all cases. Indeed, the inclusion of the transient outward current in the model allows to fit perfectly the notch in the transmembrane voltage after the initial depolarization phase.
Let us mention that human data provided by other authors[24, 27, 32, 33] differ slightly. The cause of the discrepancies being presumably due to the variation in the temperature of experiment, the heart rate stimulation, and if the right/left heart side ventricle were used. In the present paper, we considered the experiments done by Li et al.[24] because in the same paper the authors present the three different types of myocytes at several physiological heart rates.
In Figure2 we compare typical experimental conduction velocities for epicardial cells[34] and the ones obtained with our model, for two different values of the diffusion coefficient. Although the dispersion in the experimental data is significant (maybe because of different diffusivities or slightly different sodium conductance), the results provided by our model agree well with typical observed conduction velocities in epicardial cells.
Comparison with detailed ventricular electrophysiological models
Following the same protocol, we have also fitted our model to the epicardial human ventricular model developed by ten Tusscher et al.[8]. In this case, since we have access to data of both APD and CVrestitution, we perform a simultaneous optimization search of the model parameters for both curves, obtained simulating Eq. (1) in a 1D cable. Here the comparison of the AP morphology has been performed at four pacing periods, i.e., 270, 300, 500, 1000 ms. Due to the memory of the TNNP model the data of the AP comparison are obtained when the stationary state is achieved. We also include the CV at those periods in the objective function that we minimize. This is an optimization problem in a highly dimensional parameter space, with presumably many local minima. It is therefore difficult to ascertain if a given solution of the minimization procedure is the best possible one. However, as it was the case with the experimental data, the comparison between the simplified model and the TNNP model is very good, as it can be seen in Figure3. For comparison sake we also show the results obtained with the model proposed by Bueno et al.[19] (BCF model), with the parameter set they provide for their fit to the TNNP model. The large disagreement between the BCF and TNNP models suggests that there is most presumably a typo error in some of the BCF parameters[19].
Another important physiological model is the LRd model for guinea pig[4, 30] (where we consider the last version of the model, obtained from[35]). In Table1 we provide the parameter values for the fit to the LRd model. In this case, the comparison between our model and the LRd model was again done at four following pacing periods: 100, 250, 500 and 1000 ms. The LRd dynamic model (using D = 10^{−3} cm^{2}/ms) gives an approximate propagation speed of 52 cm/s at large DI. This value compares well with the actual velocity ∼50 cm/s obtained with our model with parameter values fitted to the LRd model (see Table1).
Gates dynamics and transmembrane currents
To obtain a better understanding of the newly proposed model, we will briefly discuss the shape of the different currents entering into the formulation of our model. The currents in our simplified model have all dimension of inverse of time (units ms^{−1}). In order to perform a direct comparison with the currents in more complex models they must be multiplied by the factor${C}_{m}\mathrm{\Delta}\stackrel{~}{V}$ to obtain the corresponding dimensional current in μ A cm^{−2}.

J_{ fi }: fast inward current. It corresponds mainly to the Na^{+} current in more detailed electrophysiological models. For instance, in the TNNP model[8], it is expressed by
$${I}_{\mathit{\text{Na}}}/{C}_{m}={G}_{\mathit{\text{Na}}}{m}^{3}\mathit{\text{hj}}(\stackrel{~}{V}{E}_{\mathit{\text{Na}}}),$$(8)where m is an activation gate, h is a fast inactivation gate, and j is a slow inactivation gate. The sodium conductance in the TNNP model is equal to G_{ Na }= 14. 838 nS/pF. Straightforward dimensional analysis allows to compare this value for the conductance to the one in our model g_{ fi }= 14 ms^{−1}(nS/pF ≡ ms^{−1}, see the value given in the ninth row and fourth column of Table1). For the gates dynamics, some simplifications are done in the model. The dynamics of the activation gate m is fast (in realistic models τ_{ m }∼ 0. 1 ms), so we take its steady state value, given in Figure4 (although this may lead to a decrease in the numerical accuracy of the model, see[36]). Therefore, it opens at V > V_{ c }, (for the fit to TNNP, we set V_{ c }= 0. 35, equivalent to$\stackrel{~}{V}=50$ mV). The product of the steady state of the inactivation gates h and j is modeled by a step function Θ(V−V_{ c }). The time constants for opening and closing are${\tau}_{{h}_{+}}=90.5$ ms and${\tau}_{{h}_{}}=6.61$ ms, respectively, which again are of the order of magnitude of those found in the TNNP model. The Nernst potential E_{ Na } is typically in the range E_{ Na }≃ [40−70] mV, while in our case we have set V_{ fi }= 1. 18, corresponding to$\stackrel{~}{V}=33$ mV. The time constant${\tau}_{{h}_{+}}$ is related to the recovery of the gates from inactivation, and gives the CV restitution properties of the propagating AP.

J_{ si }: slow inward current. This corresponds approximately to the Ca^{2+} current in more detailed electrophysiological models, responsible for the formation of the plateau (phase 2) in the AP. The principal component is the current I_{ CaL }, that, in the TNNP model is given by
$${I}_{\mathit{\text{CaL}}}/{C}_{m}={G}_{\mathit{\text{CaL}}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}d\phantom{\rule{0.3em}{0ex}}f\phantom{\rule{0.3em}{0ex}}{f}_{\mathit{\text{Ca}}}\frac{4\stackrel{~}{V}{F}^{2}}{\mathit{\text{RT}}}\frac{C{a}_{i}{e}^{2\stackrel{~}{V}F/\mathit{\text{RT}}}0.341C{a}_{0}}{{e}^{2\stackrel{~}{V}F/\mathit{\text{RT}}}1}$$(9)This current depends on the intra and extracellular Ca^{2+} concentrations denoted by C a_{ i }and C a_{0}, respectively. The current defined by Equation (9) possesses three gates, two voltage dependent gates d and f, and a third one f_{ Ca }that depends on the calcium concentration. In the present model we do not include the calcium concentration dynamics, and therefore all the gates are only voltage dependent. We assume that voltage activated gate d is fast, therefore we neglect its dynamics and rather take its steadystate value. Formally, it means that the dynamics of calcium current is slaved to voltage dynamics, so the last part of I_{ CaL }gives an effective voltage activated gate. In the model, the dependence on voltage is given by two sigmoidal functions that determine the voltage range at which this current is active. In conclusion, all the dynamics associated with calcium lies in the voltage inactivation gate f, that operates at a time scale of${\tau}_{{f}_{+}}\sim 43$ ms for opening, and${\tau}_{{f}_{}}\sim 181$ ms, for closing. The time constant${\tau}_{{f}_{}}$ is related with the AP duration, while${\tau}_{{f}_{+}}$ determines the AP restitution properties.

J_{ to }: transient outward current. This current corresponds to the transient K^{+} current responsible for the notch in the action potential. In the TNNP model it is expressed by :
$${I}_{\mathit{\text{to}}}/{C}_{m}={G}_{\mathit{\text{to}}}\mathit{\text{rs}}(\stackrel{~}{V}{E}_{k}).$$(10)For the reversal potential we take V_{ to }= 0, which corresponds to$\stackrel{~}{V}=85$ mV, while typical values in the literature are E_{ K }≃ −100 mV. For the gates, once again we adopt the simplified formulation through step functions for opening (or closing) at V_{ r }= 0. 6, corresponding to a potential$\stackrel{~}{V}=25$ mV, similar to the TNNP model. In our model, the dynamics of the gate s is similar to the inactivation gate of Na^{+}, and closes in a time scale of ∼4 ms. For the activating gate r, we have a time scale for opening of ∼13 ms. A comparison of the steady state values of the gates as a function of the transmembrane potential is shown in Figure4.

J_{ so }: slow outward current. This current corresponds to the slow (mainly K^{+}) repolarization currents. Among the slow potassium currents, one can distinguish: I_{ Ks }, I_{ Kr }, I_{K 1},.. but all in all, their sum is almost constant. We therefore maintain the assumption of the FK model that this current depends on a single gate at steady state.
In Figure5 we show the AP and currents for the LRd model. We have grouped the currents into four groups corresponding to the classification in our model. It should be noted that our fast inward current J_{ fi } corresponds to the sum of the sodium current I_{ Na } plus the fast part of the calcium current I_{ CaL }, in more realistic electrophysiological models, as for example in the LRd. The current corresponding to the pumps I_{ NaCa }and I_{ NaK }has been grouped with the sum of the potassium currents, that collectively correspond to the J_{ so } of the simplified model. We should stress that these currents do not have a direct electrophysiological meaning, since we are dealing with a reduced model, but a comparison with the currents in realistic models helps to clarify the meaning of the different currents in the simplified model.
Comparison with other simplified models
In this section we compare our model with other existing simplified models proposed in the literature. Two of the simplest models of cardiac dynamics contain only two variables. These are the models by Mitchell and Schaeffer[16] and by Aliev and Panfilov[14]. In both cases the equation for the transmembrane voltage includes just two currents, one outward J_{ out } and one inward J_{ in }, such that dV/dt = J_{ in }−J_{ out }. Besides the equation for the transmembrane voltage, there is an additional equation for a gate variable that, in[16] modulates the inward current, while in[14] modulates the outward current. Despite their extreme simplicity, these models present restitution properties and morphologies that resemble those of cardiac cells (except for the missing notch in phase 1 of the AP). Interestingly enough to mention, in the model by Mitchell and Schaeffer it is even possible to calculate the APD restitution curve analytically, under certain constraints in the parameters. This, together with the conduction velocity at a given pacing rate (or the maximum CV) fixes univocally all the coefficients in the model. Once this is done, there is no more freedom to fit the shape of the CVrestitution curve. Unfortunately, for typical restitution curves found in human cardiac cell, the previous constraints are not fulfilled and, thus, the simple expressions for the APDrestitution do not hold.
The MitchellSchaeffer two variable model is a simplification of the FentonKarma model[15]. In the latter the inward current is split into a fast, J_{ fi } and a slow J_{ si }component. This model is able to fit both APD and CVrestitution curves at the same time, however it does not fit properly the morphology. A further simplified model was proposed by Bueno et al[19], including an additional gate that modulates the slow inward current, mimicking the effect of the transient outward current I_{ to }. This modification implies the need to introduce different values of the parameters at different levels of the transmembrane voltage, through the inclusion of Heaviside functions. In the present paper, we follow a similar approach but, instead of modulating the slow inward current as it is done in the BCF model, we explicitly introduce a fast outward current similar to the physiological I_{ to }, an approach that we reckon more natural. This approach, for instance, has been useful to study Brugada syndrome[25, 29], where the loss of the dome is achieved just changing the conductance of this new current J_{ to }.
In Figure6 we show the APD and CV restitution curves obtained for the different simplified models. For our model and BCF we consider the values fitted to human ventricular epicardial cells, while for the models by AlievPanfilov and MitchellSchaeffer we consider typical values of the parameters cited in the original papers. The latter do not correspond to human cells.
Two dimensional simulations
In this subsection, we will check that our model reproduces the complex spiral wave dynamics observed in realistic cardiac models (for a comparison of spiral wave dynamics in several cardiac models, see[37]). In order to do so, we take a piece of uniform epicardium tissue and apply an ectopic activation to study the stability of the artificially created spiral wave and its dynamics. Spiral tips may follow different types of trajectories, from circles to flower like patterns to chaotic meanderings[38, 39]. The spiral tip is determined by the algorithm described by Fenton and Karma in the reference[15]. The wave tip is defined as the point where the excitation wavefront meets the repolarization waveback of the AP. This point (x_{ tip },y_{ tip }) is the intersection point of the lines V = V_{ iso }and ∂_{ t }V = 0. The arbitrary choice for V_{ iso }only slightly modifies the meander trajectories. In this article, we have chosen V_{ iso }as a half of the maximum depolarization potential value, i.e. V_{ iso }≈ (ΔV)/2.
In Figure7 it is shown the spiral dynamics for the LRd model and our present model. The size of the 2D domain is 10x10 cm. The time integration spans 4,000 ms and only a fraction is shown. The spirals are stable and the spiral tip experiences a meander that is comparable for both the original LRd model and our model. Here we should point out that as the spiral dynamics is not directly related to the APD and CV curves, this comparison constitutes an additional test for the validity of our model. However, given the strong memory effects in the LRd model, a very precise comparison is difficult to make, since the spiral form and dynamics is slowly changing with time. Quantitatively, the period of the spiral wave for the LRd model is T = 104. 8 ms (APD=86.8 ms and DI=18 ms) and for our model fitted for LRd, the spiral period is T = 102 ms (APD=83.8 ms and DI=18.2 ms). The high frequency of the spiral waves is the reason for using the four different BCL in the fitting of the AP morphology in section II (especially the high frequency rate at BCL=100 ms). Actually, due to the spiral motion, the period varies slightly from point to point in the two dimensional domain (Doppler shift). The size of the spiral core is also comparable and is of the order of 1cm (see Figure7).
Three dimensional simulations and ECG calculations
Cardiac arrhythmias, especially those occurring in the ventricles, are intrinsically threedimensional phenomena whereas direct experimental observations are only mostly limited to surface recordings. Other techniques like plunge electrodes can and have been used to obtain intramural information, albeit at coarse spatial resolutions[40], and transillumination is another technique that can be also used in some circumstances[41]. Therefore, it is in general difficult to identify the exact nature of the arrhythmic states and the mechanisms causing their initiations. In this respect, computer simulations of the propagating action potential throughout the heart constitute an invaluable tool to get a better insight into these mechanisms.
In this section we exploit a three dimensional computer model of the rabbit heart ventricles to further compare the original LRd model and our new model. In particular, we compute and compare the pseudoECGs[42] resulting from the AP propagation for both models. The details of the numerical scheme used in this section can be found in the paper by Bragard et al.[43]. The geometry of the ventricles are corresponds to the rabbit heart[26] and we have also included the anisotropy of the conductivity tensor which is very important in order to achieve realistic simulations of the propagating AP. Here, the numerical grid is composed by cubic voxels of 0. 025 cm size and the time step of the explicit scheme is set to dt = 0. 01 ms. Using such a small time step ensures to capture the smallest time scale linked to the depolarization of the heart myocytes (≈1 ms). The values for the longitudinal diffusion D_{∥} = 10^{−3} cm^{2}/ms and transverse diffusion D_{⊥} = 6. 75×10^{−5}cm^{2}/ms are taken from the paper by Aguel et al.[44] and reflect the fact that the conduction velocity is more than three times faster along the fibers than in the plane transverse to the fibers (i.e. transverse isotropic diffusion tensor). In addition to the geometry and the fiber anisotropy, we have modeled the complex fast conduction system (His bundle and Purkinje fibers[45]) by a time dependent external activation along the endocardial layer as it is done in the article by Boulakia et al.[46]. Indeed, the time sequence of the depolarization of the ventricular tissue is important for the reconstruction of the ECG. Here, in the rabbit heart ventricles, we have represented in black the location where the activation is initiated (see Figure8). The detail of the sequence of the firing is as follows: At t=0, the fibers between the base and the midseptum are fired; at t=5 ms, the fibers between the midseptum and the apex are fired and, finally, the remaining fibers are fired at t=10 ms. In the present paper, the ECG has been reconstructed using the heart dipole technique[47]. This technique consists in solving the forward problem of electrocardiography with a monodomain approximation[46] and assuming that the torso is an homogenous medium. This is the simplest manner to compute the pseudoecgs but still gives satisfactorily qualitative results. The method consists in adding all the microscopic dipoles (created at the depolarization fronts) into a single vector which is called the heart dipole vector as follows:
where in Equation (11),$\mathcal{V}$ represents the integration volume (ventricles) and$\overrightarrow{r}$ is the vector joining the geometrical center of the heart to the microscopic dipoles created by the moving fronts. The next step consists in projecting the heart vector on some standard directions in order to compute the different derivations of the ECG. In our case, because of the crude approximation used for the torso, we are unable to compute the chest leads[46]. As the main interest of this section was to compare our model with the LRd model, we will only compute the standard lead I and check if they are qualitatively the same.
For the comparison of the pseudoecgs we have stimulated both our model and the LRd model at a period of BCL=400 ms during several beats. In the case of the LRd model, and because this model is known to have memory effects associated to the slow dynamics of ion concentrations, we have started the simulation from an initial condition obtained from a separate one dimensional simulation of one thousand beats. The electrical organized activity associated with the AP propagation is well captured in the corresponding pseudoECGs as shown in Figure9. We observe that after few beats both models (ours and LRd) converge as expected. It should be noticed that the Twave indicating the repolarization phase of the ventricles is inverted in both ECGs in Figure9. This is due to the fact that our description of the ventricles lacks to include the transmural heterogeneities, i.e. different APDs of the cardiac myocytes in the epi, mid and endocardium tissues[46, 48]. From the computational speed point of view, the difference between our model and the LRd model is striking. In order to have a fair comparison, we have simulated both models on the same single processor, thus avoiding differences linked to parallelization speed variations. For the LRd model the simulation lasted for 68 hours, while only 4h30’ for our model. Therefore a gain of approximately a factor twenty in computational speed. As mentioned in[46], if we want to address the inverse problem of electrocardiography the speed is a crucial factor to take into account because one has to solve many times the forward problem to get an approximation of the inverse problem. The use of a simplified model is one way to deal with the speed issue.
Conclusions and further works
In this paper we have analyzed a semiphysiological cardiac model that has been used previously to study reexcitation in a medium with dispersion of repolarization[25]. The model reproduces well the action potential morphology of different cardiac cells from experimental data and also from more complete physiological models.
The CVrestitution and the pseudoECG’s in a realistic model of heart ventricles have also been satisfactorily compared. One concludes that to study AP propagation our present simplified model is a good alternative to the most costly physiological models. The gain in computer speed is around a factor of twenty and the need for computer memory (RAM) is also greatly reduced. In certain situations it also allows to obtain analytically some characteristics of propagation[25], which, in a realistic cardiac model, would be intractable. A limitation of the model is the absence of memory due to the lack of the slow dynamics of ion concentrations. This makes unfeasible to fit simultaneously dynamic and S1S2 restitution curves and may not give a good comparison with realistic models under situations where the stimulation frequency is abruptly changed. In the near future, we plan to refine the present three dimensional simulation to include the important transmural heterogeneities[49], in order to study transmural reexcitations and to obtain more realistic ECGs from the 3D model.
References
 1.
Hodgkin A, Huxley A: A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol. 1952, 117: 500544.
 2.
Fenton FH, Cherry EM: Models of cardiac cell. Scholarpedia. 2008, 3: 186810.4249/scholarpedia.1868.
 3.
Puglisi JL, Bers DM: LabHEART: an interactive computer model of rabbit ventricular myocytes ion channels and Ca transport. Am J Physiol. 2001, 281: C2049C2060.
 4.
Luo C, Rudy Y: A dynamic model of the cardiac ventricular action potential  Simulations of ionic currents and concentration changes. Circ Res. 1994, 74: 10711097. 10.1161/01.RES.74.6.1071.
 5.
Pandit SV, Clark RB, Giles WR, Demir SS: A mathematical model of action potential heterogeneity in adult rat left ventricular myocytes. Biophys J. 2001, 81: 30293051. 10.1016/S00063495(01)759437.
 6.
Winslow RL, Rice J, Jafri S, Marban E, O’Rourke B: Mechanisms of altered excitationcontraction coupling in canine tachycardiainduced heart failure. II. Model studies. Circ Res. 1999, 84: 571586. 10.1161/01.RES.84.5.571.
 7.
Priebe L, Beuckelmann DJ: Simulation study of cellular electric properties in heart failure. Circ Res. 1998, 82: 12061223. 10.1161/01.RES.82.11.1206.
 8.
ten Tusscher K, Noble D, Noble PJ, Panfilov AV: A model for human ventricular tissue. Am J Physiol Heart Circ Physiol. 2004, 286: H1573—H1589
 9.
Iyer V, Mazhari R, Winslow RL: A computational model of the human leftventricular epicardial myocyte. Biophys J. 2004, 87: 15071525. 10.1529/biophysj.104.043299.
 10.
Grandi E, Pasqualini FS, Bers DM: A novel computational model of the human ventricular action potential and Ca transient. J Mol Cell Cardiol. 2010, 48: 112121. 10.1016/j.yjmcc.2009.09.019.
 11.
Clancy CE, Rudy Y: Na+ channel mutation that causes both Brugada and longQT syndrome phenotypes: a simulation study of mechanism. Circulation. 2002, 105: 12081213. 10.1161/hc1002.105183.
 12.
FitzHugh R: Impulses and physiological states in theoretical models of nerve membrane. Biophysical J. 1961, 1: 445466. 10.1016/S00063495(61)869026.
 13.
Nagumo J, Arimoto S, Yoshizawa S: An active pulse transmission line simulating nerve axon. Proc IRE. 1962, 50: 20612070.
 14.
Aliev RR, Panfilov AV: A simple twovariable model of cardiac excitation. Chaos: Solitons Fractals. 1996, 7: 29310.1016/09600779(95)000895.
 15.
Fenton F, Karma A: Vortex dynamics in three dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation. Chaos. 1998, 8: 2047. 10.1063/1.166311.
 16.
Mitchell CC, Schaeffer DG: A twocurrent model for the dynamics of cardiac membrane. Bull Math Biol. 2003, 65: 767793. 10.1016/S00928240(03)000417.
 17.
Cherry EM, Fenton FH: Suppression of alternans and conduction blocks despite steep APD restitution: electrotonic, memory, and conduction velocity restitution effects. Am J Physiol Heart Circ Physiol. 2004, 286: H2332—H2341
 18.
Cherry EM, Ehrlich JR, Nattel S, Fenton FH: Pulmonary vein reentry–properties and size matter: insights from a computational analysis. Heart Rhythm. 2007, 4: 15531562. 10.1016/j.hrthm.2007.08.017.
 19.
BuenoOrovio A, Cherry EM, Fenton FH: Minimal model for human ventricular action potential in tissue. J Theor Biol. 2008, 253: 544560. 10.1016/j.jtbi.2008.03.029.
 20.
Cytrynbaum E, Keener JP: Stability conditions for traveling pulse: modifying the restitution hypothesis. Chaos. 2002, 12: 788799. 10.1063/1.1503941.
 21.
Antzelevitch C: Role of spatial dispersion of repolarization in inherited and acquired sudden cardiac death syndromes. Am J Physiol Heart Circ Physiol. 2007, 293: H2024H2038. 10.1152/ajpheart.00355.2007.
 22.
Brugada P, Brugada J: Role of spatial dispersion of repolarization in inherited and acquired sudden cardiac death syndromes. J Am Coll Cardiol. 1992, 20: 13911396. 10.1016/07351097(92)90253J.
 23.
Jervell A, LangeNielsen F: Congenital deafmutism, functional heart disease with prolongation of the QT interval and sudden death. A Heart J. 1957, 54: 5968. 10.1016/00028703(57)900790.
 24.
Li GR, Feng J, Yue L, Carrier M: Transmural heterogeneity of action potentials and Ito1 in myocytes isolated from the human right ventricle. Am J Physiol. 1998, 275: H369H377.
 25.
Cantalapiedra IR, Peñaranda A, Echebarria B, Bragard J: Phase2 reentry in cardiac tissue: role of the slow calcium pulse. Phys Rev E. 2010, 82: 011907
 26.
Vetter FJ, McCullogh AD: Threedimensionnal analysis of regional cardiac function: a model of the rabbit ventricular anatomy. Prog Biophys Mol Biol. 1998, 69: 157183. 10.1016/S00796107(98)000066.
 27.
Nabauer M, Beuckelmann DJ, Uberfuhr P, Steinbeck G: Regional differences in current density and ratedependent properties of the transient ventricular electrophysiology. Am J Physiol. 1996, 292: H43—H55
 28.
Cantalapiedra IR, Peñaranda A, Mont L, Brugada J, Echebarria B: Reexcitation mechanisms in epicardial tissue: role of I(to) density heterogeneities and I(Na) inactivation kinetics. J Theor Biol. 2009, 259: 850859. 10.1016/j.jtbi.2009.04.021.
 29.
Peñaranda A, Cantalapiedra IR, Echebarria B: Slow pulse due to calcium current induces phase2 reentry in heterogeneous tissue. Comput Cardiol. 2010, 37: 661664.
 30.
Faber GM, Rudy Y: Action potential and contractility changes in Na+ overloaded cardiac myocytes: a simulation study. Biophys J. 2000, 78: 23922404. 10.1016/S00063495(00)76783X.
 31.
Webpage with the codes used for the fit and further information.http://wwwfa.upc.es/websfa/eupb/NOLIN/CARDIAC/Simp_model.html,
 32.
Drouin E, Charpentier F, Gauthier C, Laurent K, LeMarec H: Electro physiologic characteristics of cells spanning the left ventricular wall of human heart: evidence for presence of M cells. J Am Coll Cardiol. 1995, 26: 185192. 10.1016/07351097(95)00167X.
 33.
Morgan JM, Cunningham D, Rowland E: Dispersion of monophasic action potential duration: demonstrable in humans after premature ventricular extrastimulation but not in steady state. J Am Coll Cardiol. 1992, 19: 12441253. 10.1016/07351097(92)90331G.
 34.
Yue AM, Franz MR, Roberts PR, Morgan JM: Global endocardial electrical restitution in human right and left ventricles determined by noncontact mapping. J Am Coll Cardiol. 46: 10671075.
 35.
LRd ventricular cell model (guineapigtype), source code. [http://rudylab.wustl.edu/research/cell/methodology/cellmodels/LRd/code.htm], []
 36.
Bernus O, Wilders R, Zemlin CW, Verschelde H, Panfilov AV: A computationally efficient electrophysiological model of human ventricular cells. Am J Physiol Heart Circ Physiol. 2002, 282: H2296—2308
 37.
Shajahan TK, Nayak AR, Pandit R: Spiralwave turbulence and its control in the presence of inhomogeneities in four mathematical models of cardiac tissue. PLoS ONE. 2009, 4: e473810.1371/journal.pone.0004738.
 38.
Davidenko JM, Pertsov AM, Salomonsz R, Baxter WT, Jalife J: Stationary and drifting spiral waves of excitation in isolated cardiac muscle. Nature. 1992, 355: 349351. 10.1038/355349a0.
 39.
Fenton FH, Cherry EM, Hastings HM, Evans SJ: Multiple mechanisms of spiral wave breakup in a model of cardiac electrical activity. Chaos. 2002, 12 (3): 852892. 10.1063/1.1504242.
 40.
Allison JS, Qin H, Dosdall DJ, Huang J, Newton JC, Allred JD, Smith WM, Ideker RE: The transmural activation sequence in porcine and canine left ventricle is markedly different during longduration ventricular fibrillation. J Cardiovasc Electrophysiol. 2007, 18: 306312.
 41.
Baxter WT, Mironov SF, Zaitsev AV, Jalife J, Pertsov AM: excitation waves inside cardiac muscle using transillumination. Biophys J. 2001, 80: 516530. 10.1016/S00063495(01)760341.
 42.
Bernabeu MO, Corrias A, PittFrancis J, Rodriguez B, Bethwaite B, Enticott C, Garic S, Peachey T, Tan J, Abramson D, Gavaghan D: Grid computing simulations of ion channel block effects on the ECG using 3D anatomicallybased models. Comput Cardiol. 2009, 36: 213216.
 43.
Bragard J, Marin S, Cherry E, Fenton F: Validation of a model of cardiac defibrillation. To appear in Springerbook (2012)
 44.
Aguel F, Eason J, Trayanova N: Advances in modeling cardiac defibrillation. Int J Bifurcation Chaos. 2003, 13: 37913805. 10.1142/S0218127403008892.
 45.
Katz AM: Physiology of the heart. 2005, Philadelphia: Lippincott Williams & Wilkins
 46.
Numerical simulations of electrocardiograms. Edited by: Ambrosi D, Quarteroni A, Rozza G. 2011, Springer
 47.
Constanzo LS: Physiology. 2002, Philadelphia, P A: W.B. Saunders
 48.
Gussak I, Antzelevitch C: Cardiac Repolarization: Bridging Basic and Clin Sci. 2003, New York: Humana Press, Springer
 49.
Abd Allaha ESH, Aslanidic OV, Telleza JO, Yannia J, Billeterd R, Zhangc H, Dobrzynskia H, Boyett MR: Postnatal development of transmural gradients in expression of ion channels and Ca2+ handling proteins in the ventricle. J Mol Cell Cardiol. 2012, 53 (2): 145155. 10.1016/j.yjmcc.2012.04.004.
Acknowledgements
This work has been financially supported by MICINN (Spanish Ministry of Science) under coordinated projects FIS200806335 and FIS201128820. The work by AP has also been supported under project FIS2012376550201. The authors acknowledge useful discussions with Professors S. Santidrian and I. GarciaBolao (UNAV) about some physiological questions. Discussions about the model and its simulation with Drs. F. Fenton, E. Cherry and N. Otani are also acknowledged.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
All authors read and approved the final manuscript.
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/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Received
Accepted
Published
DOI
Keywords
 Action Potential Duration
 Spiral Wave
 Action Potential Propagation
 Pace Rate
 Transmembrane Voltage