Mathematical Problems in Engineering
Volume 2020 (2020), Article ID 3704523, 12 pages
https://doi.org/10.1155/2020/3704523

Research Article

A Graphical-Based Mathematical Model for Simulating Excitability in a Single Cardiac Cell

S. H. Sabzpoushan and A. Ghajarjazy

Department of Biomedical Engineering, Iran University of Science and Technology (IUST), Tehran 16846-13114, Iran

Correspondence should be addressed to S. H. Sabzpoushan; sabzposh@iust.ac.ir

Received 27 December 2019; Revised 1 June 2020; Accepted 12 June 2020; Published 22 July 2020

Academic Editor: Alessio Gizzi

Copyright © 2020 S. H. Sabzpoushan and A. Ghajarjazy. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Excitability is a phenomenon seen in different kinds of systems, e.g., biological systems. Cardiac cells and neurons are well-known examples of excitable biological systems. Excitability as a crucial property should be involved in mathematical models of cardiac cells, along with the other biological properties. Excitability of mathematical cardiac-cell models is usually investigated in the phase plane (or the phase space) which is not applicable with simple mathematical analysis. Besides, the possible roles of each model parameter in the excitability property cannot be investigated explicitly and independently using phase plane analysis. In this paper, we present a new graphical-based method for designing excitability of a single cardiac cell. Each parameter in the presented approach not only has electrophysiological interpretation but also its role in regulating excitability is evident and can be analysed explicitly. Our approach is simpler and more tractable by mathematical analysis than the phase plane method. Another advantage of our approach is that the other important feature of the cardiac cell action potential, i.e., plateau morphology, can be designed and regulated separately from the excitability property. To evaluate our presented approach, we applied it for simulating excitability in well-known complex electrophysiological models of ventricular and atrial cells. Results show that our model can simulate excitability and time evolution of the plateau phase simultaneously.

1. Introduction

Mathematical modelling of cardiac myocytes proves to be a useful tool for investigating the dynamics of electrical excitation and contraction in the heart [1, 2]. The mathematical description of cells was initiated by the pioneering work of Hodgkin and Huxley (HH model) [3], who first described quantitative analysis of the ionic currents of the squid giant axon. Ten years later, Noble developed a mathematical model based on the HH model to gain insights into the dynamics of the Purkinje cell action potential (AP) [4].

Later on, ionic models that describe the cardiac action potential have become increasingly more complex and more realistic [5–7]. However, the complexity of these models makes it difficult to identify a subset of responsible parameters for an observed behaviour. Since we are looking for approaches which are tractable by simple mathematical analysis, one way for avoiding complex models is to present simple procedures and construct simplified models of excitable media. As an example of simple approaches and models, we can name the FitzHugh and Nagumo (FHN) model [8–10]. FHN provides a qualitative description of generic excitable media, as well as information on effects of a reduced set of parameters. Fenton–Karma (FK) model [10] is a widely used simplified model for cardiac AP, which has three state variables. The parameters of the FK model can be chosen to reproduce some properties of ionic models and experimental data for different cardiac cells [11–13].

In electrophysiological models, the parameters often have physiological meaning and therefore are typically obtained by fitting the model to experimental data [14]. Many studies have been published on different methods for tuning models to fit desired features of single AP [15–19]. However, the fitted parameter set will almost certainly not be the best possible set [20–22].

Excitability is a basic well-known phenomenon seen in biological systems, e.g., cells and neurons. Cardiac myocytes consist of excitable cells that can generate and propagate excitations. Excitability in cells is due to the electrophysiological behaviours of ion channels and currents [23].

Figure 1(a) shows a real AP of a ventricular cell (VTC) [5], where the dark line represents the AP. In this figure, we have segmented the dark line into two types of dashed and dotted lines. In dashed-line sections, we see upstrokes (depolarization), downstrokes (repolarization), and silence (rest). A look at the electrophysiological literature shows that the dashed-line segments of the AP, shown in Figure 1(a), reflect the excitability of a cell [23–28], which is the purpose of the study.

(a)

(b)

(c)

(d)

Figure 1: Schematic review of our idea and main definitions in this research. (a) How a real ventricular cell action potential can be assigned to our proposed formalism of excitability. (b) Definition of excitability features. (c) Phase plane representation of an excitable system. Phase plane analysis has limitations and burdens. (d) General formalism for excitable system modelling (see text for details).

In this research, an excitable cell (EXC) is defined as a system that has only one rest state, and if it perturbs from the rest, according to small (subthreshold) or large (superthreshold) perturbations, it will take two different trajectories for returning to the rest state. In other words, this system has a single stable equilibrium point but two modes of returning to that point [23, 24]. In this paper, we are interested in cell membrane voltage, , as the concerned state variable of the EXC. Figure 1(b) illustrates excitability features of an EXC generally, i.e., definitions and quantities, which introduce a given VTC as an EXC. Regarding Figure 1(b), these features are , ,, , , and TE. Here, stands for the rest state or rest voltage, is the upstroke threshold, i.e., if the applied stimulus to the EXC puts the membrane voltage below , the EXC will not be excited, and the membrane voltage, , will return to , and if the applied stimulus locates the membrane voltage above it, the EXC will be excited, and the membrane voltage experiences a large evolution, going to the peak state, ; and then, after a time evolution TE (with the duration of ), it returns to the rest state, , through . Here, is a point at which fast downgoing (downstroke) of membrane voltage initiates.

The main aim of this paper is to present a mathematical method for modelling the excitability property of a single cardiac VTC. Our approach in the presented method is graphical, besides each parameter of the model not only has an electrophysiological interpretation but also its effect on excitability can be analysed independently.

2. Materials and Methods

In different panels of Figure 1, we have summarized schematically the main points and definitions of our method and model. Figure 1(a) illustrates how our mathematical framework in Figure 1(b) may be adapted to a real ventricular action potential (VAP). It is seen that a VAP is composed of five phases [29, 30]. Phase 4 represents the rest state of the VAP, where the value of the cell membrane potential is stable, . Phase 0 or depolarization phase is the counterpart of the fast upstroke phase in Figure 1(b), in which excitability threshold voltage () should be determined experimentally. In phase 1, a spike is generated in the VAP so that the spike voltage can be treated as in Figure 1(b). Phase 2 or plateau phase is an important part in which the VAP usually experiences a slow repolarizing time evolution. It is evident that we may regard phase 2 as TE in Figure 1(b). Phase 3, which is called the repolarization phase formally, is in fact a fast repolarization or downstroke of the AP. As it is seen in Figure 1(a), may be regarded as an edge point between phases 2 and 3 at which a slow repolarization switches to a fast one. Mathematically, this point may be defined as the right-hand intersection point of dashed and dotted lines in Figure 1(a), where the slope of the right dashed line is equal to the maximum slop of repolarization, and the dotted line is a line drawn from the tangent to the plateau. In Figure 1(a), action potential duration (APD) is the time interval from the onset of depolarization to the point of 40%, 50%, or 90% repolarization. This characteristic is denoted as APDn (n = 40, 50, and 90) [31]. In Figure 1(b), it is plausible to assume as APD in the VAP.

Phase plane (space) approach is a common method for analyzing dynamics of neurons and other excitable cells [26]. Although excitability is one of the dynamical manifestations of a cell, phase space analysis suffers from some limitations and burdens in this context. To be more clear, an abstract definition of excitability dynamics of a cell in the phase plane is illustrated in Figure 1(c) [26]. In Figure 1(c), a long trajectory starts near a small square and then returns to the small dark point. Indeed, the large trajectory is an AP, and the small dark point is the rest state. Every trajectory that starts in and around the small neighborhood of the square leaves it and then returns to the dark point (rest). The vector fields in the small dashed circle as well as outside of it have important impacts not only on excitability but also on morphology and other properties of the generated AP. These vectors are affected by almost all system parameters which make mathematical analysis of such a system difficult in the phase plane. Especially, the effects of parameters on excitability and AP morphology cannot be analysed individually. Moreover, parameters responsible for depolarization, plateau, and repolarization phases cannot be adjusted separately. In Figure 1(c), large trajectory intersects itself, so it may need to analyse the system in more than two state spaces, which make the mathematical analyzing more complex.

As we will see, our graphical-based mathematical approach is not only simple and analyzable but also the role of each parameter in system properties can be investigated directly and almost independently. Moreover, analyzing of systems with more than two state variables is expected to be simple and straightforward.

From electrophysiology point of view, a cell can be regarded as a capacitor which is charged (or discharged) by different ionic currents, where ionic currents are generated by a variety of mechanisms [26]. Figure 1(d) illustrates schematically the above fact, where k is the number of currents. This figure shows the main formalism of all Hodgkin–Huxley- (HH-) based VTC models [3], as well as our model formalism.

2.1. Graphical Approach and Model Equations

Regarding Figure 1(d) and general formalism of HH, our observation from a VTC is a voltage which is described by the following differential equation:

In this research, we look at the above currents from a novel point of view to categorize them into two electrophysiological groups: 1-currents which may be treated as functions of the membrane voltage, with the symbolic notation for this group and 2-those which may be treated as functions of voltage and time simultaneously, with the symbolic notation . In our simple model, we use one delegate from each group with symbols and , respectively. Also, we refer to them as and for simplicity. So, we can write

where stands for the membrane capacitance and represent two types of ionic currents, sometimes named as fast and slow, respectively. The explicit expressions for these currents are given by

In equation (3), is a voltage-dependent membrane conductance, and denotes a constant voltage (Nernst voltage by some ways). In equation (4), is a voltage- and time-dependent (dynamic) membrane conductance, and is another constant voltage (Nernst voltage). In this respect, all elements of our model have exact electrophysiological interpretations. Moreover, to be more electrophysiological realistic, we write

where is a constant conductance and delegates a voltage-dependent dynamic gating variable. The following equation describes the nonlinear differential equation of this gating variable:

where , the steady state value of at voltage , is considered as a sigmoid function of as follows:

In order to prepare more adaptability of our model with electrophysiology of VTCs, we decompose as in which is a constant conductance, and is a voltage-dependent gating variable. We also consider a sigmoid function of for as

Equations (2)–(8) represent the essential formalism and structure of our two state-variable electrophysiological VTC model [1, 3–5, 10, 12, 13].

2.1.1. Graphical Approach

In a simple model which is described by equations (2)–(8), ten parameters, , should be estimated so that they can simulate a real AP of a single VTC. It is desired that the last three parameters, , which are related to TE or in some ways, have no (or less) contributions to excitability reproduction of the model.

Our goal is to tune the model (equations (2)–(8)) with the use of a graphical-based algorithm so that desired excitability properties, , are fulfilled, and the effect of each parameter on these properties can be analysed explicitly.

To present our novel graphical approach, we plot and versus on the same -axis in Figure 2. Regarding equation (2), their intersection points mean the current balance or equilibrium points.

Figure 2: A typical diagram of currents and versus . On , two graphs have intersections at three points: A, B, and D. At , the line of rotates clockwise around pivot point A. Therefore, B and D move toward each other and meet at , where the line is tangent to the graph.

Current can be drawn easily as a function of , i.e., versus . According to equation (3), it intersects the -axis at two points: and where . Interestingly, the current at each time () can be supposed as a line with the slope that intercepts-axis at . As we mentioned above, the possible intersection points of the two graphs at each indicate the equilibrium points of the system, i.e., .

Regarding the above discussion, as time increases, if be about increasing, the line rotates clockwise around its pivot point, , and if be decreasing, the line rotates counterclockwise. We see that the pivot point (that we will set in future) is a key point.

Now, suppose we have adjusted the model parameters at so that and have three intersection points (A, B, and D), where , , and as depicted in Figure 2. All of these points are equilibrium points because the equality condition is satisfied () at all of them. Assume that the model is at the rest state, i.e., ; we apply a subthreshold excitation that increases just a little () so that in this situation. It is clear from Figure 2 that ; consequently, regarding equation (2), the membrane potential will be pushed back to the rest state, . All of these mean that the model will not be excited by a subthreshold excitation.

Differently, if the model is excited by a superthreshold stimulus, then increases to so that . Regarding Figure 2 and equation (2), leads to the increase in voltage. As a consequence, the membrane potential will be attracted by the next equilibrium point . All of these mean that the model will be excited by a superthreshold excitation.

The above paragraphs reveal that the graphical illustration in Figure 2 can satisfactorily demonstrate the model of equations (2)–(8). In this approach, the model dynamics can be understood easily by a visual illustration.

(1) Parameter Estimation Algorithm. To adjust our model (equations (2)–(8)) to produce a desired AP, we should solve an inverse problem: given characteristics of an AP as were introduced in Figures 1(a) and 1(b), compute model parameters. In this way, our approach will be graphical-based, i.e., we use Figure 2 framework. As it can be deduced from Figures 1(a)–1(c), the dynamics of the AP of a VTC can be studied in two parts: excitability and TE. Since excitability is the main subject of this paper, we aim to represent the algorithm so that one can analyse these two parts as independently as possible. In the following, we show how our algorithm and graphical approach can compute model parameters in excitability and TE parts independently.

(a)Parameter computation in the excitability section: in this section, membrane voltage grows from to , through , and develops from to . A visual deduction from Figure 2 is

 Since time has not an important role in the excitability section, we can consider that all things happen instantaneously. Therefore, has not been activated throughout this section and remains at its value before the stimulated state, . Assuming for simplicity, we rewrite equations (2) and (6), respectively, as

 Current in equation (10) has two important roles. First, it makes a positive instantaneous feedback and a fast regenerative process for pushing from to . Second, as we will see later, it opposes in the plateau phase ( evolution from toward ) and controls the decrements of as well. Concerning the excitability definition, the amplitude of should be small on the to interval to suppress the subthreshold stimulus and large on the to interval in order to cause fast regenerative process. The above required properties for are provided by the definition of in equation (8). Assuming , we can write equation (12) in the triangle with corners , and :

 We define parameter α as the ratio of m-gate in two distinct voltages ( and ) as follows:

 Substituting equation (13) into equation (3), one gets

 Using equation (14), we can write

 Solving algebraic equations (12) and (15) for , one gets

 (16) means that, by having desired characteristics, , and , for a cell and choosing , we can get one of the model parameters, . Considering line, at the fixed point, we have

 For simplicity and without loss of generality, we can scale currents so that

 By substituting equations (4) and (18) into equation (17), we get

 Therefore,

 Writing current balance at another equilibrium point () leads to

 Using equations (15), (18), and (20) in equation (21), we get

 Solving equation (22) for gives

 After some algebraic manipulations on equation (23), we get equation (24), which is a counterpart of equation (16):

 Satisfaction of equations (16) and (24) simultaneously results in the following equation:

 It means that another parameter of the model, , is determined. Writing current balance at stable fixed point leads to ; therefore, we can treat α as a free parameter and choose it small enough to satisfy . According to the -gate definition, if , then which leads to , so the hope is that the rough condition be satisfied better. In this condition,

 By some algebraic manipulations on equation (26), we get a relation between parameters of -gate:

 In addition, the instantaneous process on the to interval can be convincing by defining as a saddle fixed point at time . It means that should be tangent to at point , which is fulfilled by the following condition:

 By substituting equation (27) into equation (28), we get

 Solving equation (29) leads to the determination of , which is one of the model parameters. Since it seems difficult to find an explicit solution to equation (29), we solve it graphically by drawing both the left- and the right-hand side of equation (29) on the same coordinates versus , and the cross section point is the solution. Afterward, using obtained in equation (27) results in the determination of . Now, having and , values of -gate can be found at any voltage. Therefore, can be determined by the following equation:

 In conclusion, we see that the parameters , , , , and can be computed in the excitability section with the use of features of the AP of given VTC.(b)Parameter computation in the TE section: according to our proposed approach and formalism, the TE of an AP is modelled by rotation of about its pivot point . This rotation is attributed to gradual activation of , which causes a time course of voltage progress from to . It is noteworthy that although is considered as a single state variable, it may delegate more than one state variable actually (we have deliberated this matter more in the discussion section). Generally, TE may be treated as a repolarization process as it is seen in Figure 1(a).

Here, just as an example, we show how one can use our graphical approach for simulating the TE section. We assume as an activation gate and choose an arbitrary small value for it at the rest state:

Using equation (31), satisfying in equation (20), we can determine :

Recalling our graphical approach, the line should be tangent to the graph at , so

Then, by considering equations (3) and (4), one gets

Having two values of , in equations (31) and (34), we can determine and as follows:

Having in mind the above discussions about excitability and TE sections and considering Figure 2, as well as presented equations, we implemented numeric calculations of our graphical approach in Matlab as follows.

At the first time step after excitation, raises toward ; this increment makes a clockwise rotation of the line, i.e., generates and forms a new fixed point (e.g., point d shown in Figure 2). This new fixed point creates new , as well. As time progresses, a sequence of and consequently are generated. To satisfy TE of wanted AP and inducing fast repolarization at , the sequence of the generated two fixed points (d, b) should meet each other at time and at voltage . When this happens, the regenerative process occurs to get the voltage back to quickly. Regarding equation (6), the time step should be small enough so that we can capture the dynamics of during the rotation of the line. In Figure 3, the general algorithm for the determination of the successive intersection points of the line with the graph is summarized.

Figure 3: The algorithm for the determination of successive intersection points of the line with the graph.

In Figure 4, the whole steps and sequences of parameter estimation of our approach are illustrated.

Figure 4: The algorithm for implementation of our graphical approach. The desired excitability features are inputs to the algorithm, and the model parameters are its outputs.

3. Results and Discussion

In this paper, we considered the dynamics of the AP of a VTC in two divisions: excitability and TE. We presented a model and showed how its parameters can be estimated by using a graphical approach. Our method is graphical, mathematically tractable, and mostly based on algebra which is simpler than differential equations and phase plane analysis approaches [18]. Here, we remind that our main goal was to show how our graphical approach can simulate the excitability property of an AP. To examine this ability empirically, we used two complex electrophysiological models by employing OpenCOR [32], which is an electrophysiological modelling medium, and CellML repository [33]. Based on the illustrated algorithm in Figure 4, we developed a computerized algorithm in Matlab (MathWorks, version 7.10, 2010) for implementing our approach.

Table 1 lists the computed parameter values of the proposed model for reproducing excitability characteristics of Noble-2001 VTC [34], and Figure 5(a) summarizes the graphical approach for this case. The successive cross marks in Figure 5(a) indicate sequential cross sections of the line with the curve, and the cross sections at 50, −28, −49, and −92.8 mv voltages correspond to , , and , respectively. In Figure 5(b), the simulated AP is compared with the Noble-2001 model. In Table 2, we have compared the features of the examined and simulated models quantitatively.

Table 1: The estimated parameter values of the proposed model for simulating excitability features of the Noble-2001 model.

Parameters

0.0209
0.007
0.07
97.89
−92.8
4.96
−38.1
28.8
−29.43
75.5

(a)

(b)

Figure 5: Reproducing AP excitability generated by the Noble-2001 model. (a) Proposed graphical approach. (b) Simulated AP.

Table 2: Comparing the excitability features of the AP simulated by the proposed model with the ones generated by the Noble-2001 model.

T

Noble-2001

−92.8 mv
−49 mv
−28 mv
50 mv
160 msec
—
—

Parameters

−92.7
−49
−28
49
160
75.43
−0.000064

−92.55
−48.91
−28
49
160
80.52
−0.000349

In another trial, we examined the CRN-1998 atrial cell model [35]. Estimated value of parameters in our model is listed in Table 3, and graphical results are depicted in Figure 6. Table 4 compares the excitability features of the AP of CRN-1998 with those of our proposed approach; one can see that our approach is successful to reproduce the excitability of an atrial cell as well as a ventricular one.

Table 3: The estimated parameter values of the proposed model for simulating excitability features of the CRN model.

Parameters

0.0160
0.0095
0.1
87.0659
−81.19
4.1416
−47.7800
21.3548
−34.2686
135

(a)

(b)

Figure 6: Reproducing excitability of the AP generated by the CRN-1998 model. (a) Proposed graphical approach. (b) Simulated AP.

Table 4: Comparing the excitability features of the AP simulated by the proposed model with the ones generated by the CRN model.

T

CRN

−81.19 mv
−56.88 mv
−40 mv
24.49 mv
215.6 msec

Parameters

−81.19
−56
−40
23.8
214.3

Reproduced features of the above two investigated examples demonstrate the empirical validation of our proposed graphical approach. Although the main aim of this paper was modelling the excitability of the VTC using the proposed approach, as we saw in the preceding example, it can also simulate excitability in atrial cells; in fact, our approach can be applied to every excitable system which can be casted to Figure 1 formalism.

In our paper, the focus of attention was on excitability simulation, i.e., excitability section of Figures 1(a) and 1(b), and the work on the TE issue can be the matter of an independent research. However, we present a brief argument about this topic here as a preface to future research studies.

In Figures 5(b) and 6(b), although the excitability features have been provided, TEs do not seem satisfactory. .

In the above examinations, we have considered first-order dynamics for , but our approach is flexible so that we can consider different speed profiles for . To be more clear, we may even attribute higher-order dynamics to in Figure 2, i.e., rotation of about its pivot point may have higher-order dynamics. In Figures 7 and 8, we have clarified this idea by considering a more complex speed profile for the rotation of the line. In Figure 7(a), a TE morphology of the VTC is shown that is consistent with the one generated by a complex electrophysiological model [5], where in Figures 7(b)–7(d), the general considered dynamics for are illustrated, respectively. In Figure 8(a), we have done the same for Noble-2001 [34], where all values of parameters are the same as Table 1, except considering second-order dynamics for as depicted in Figures 8(b)–8(d).

(a)

(b)

(c)

(d)

Figure 7: Using higher-order dynamics for to generate a more complex morphology of the AP. (a) Generation of an AP similar to a complex electrophysiological model [5]; dashed and bold lines represent complex and simplified models, respectively. (b) time course profile. (c) speed profile. (d) acceleration profile.

(a)

(b)

(c)

(d)

Figure 8: Using sample second-order dynamics for to generate morphology of the AP of Noble-2001. (a) Generated AP; dashed and bold lines correspond to Noble-2001 and simplified models, respectively. (b) time course profile. (c) speed profile. (d) acceleration profile. The parameter values of the excitability section are the same as Table 1.

Preserving our presented method as the main framework, we can extend it to models with more than two state variables. To put it another way, we see terms like , j, , … in the electrophysiology literature [3, 5, 12, 26].In these cases, we can assume a term like as a single compound gating variable with complex dynamics, i.e., with time course like in Figure 7(b).

As a limitation of our work, we have not involved the recovery phase of in our approach still, where it has impacts on other properties of the AP like alternans [28]. In our approach in Figure 2, recovery can be defined by the mechanism of the rotation of the line after leaving tangency at the point, and and its counter clockwise rotation toward D point. This is a remaining issue which will be presented in future papers.

4. Conclusions

Mathematical modelling of ventricular cells is of great help to cardiac studies. Researchers are looking for mathematical approaches and models which are simple and tractable by mathematical analysis. In this research, we proposed a novel graphical approach which can model excitability and time evolution of the AP of a VTC independently. This approach is tractable by simple mathematical analysis, and the effect of each parameter on excitability and TE can be analysed and investigated separately. Although the approach was examined in a two state-variable model, it may be applicable for higher-order models as it was shown in the previous section.

Data Availability

No data were used to support this study.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Supplementary Materials

Visual presentation of our graphical approach, Figure 7(a).In this animation, the time evolution of currents and their intersection points are displayed. (Supplementary Materials)

References

J. L. Puglisi, F. Wang, and D. M. Bers, “Modeling the isolated cardiac myocyte,” Progress in Biophysics Molecular Biology, vol. 85, no. 2-3, pp. 163–178, 2004.

Y. Rudy and J. R. Silva, “Computational biology in the study of cardiac ion channels and cell electrophysiology,” Quarterly Reviews of Biophysics, vol. 39, no. 1, pp. 57–116, 2006.

A. L. Hodgkin and A. F. Huxley, “A quantitative description of membrane current and its application to conduction and excitation in nerve,” The Journal of Physiology, vol. 117, no. 4, pp. 500–544, 1952.

D. Noble, “A modification of the Hodgkin-Huxley equations applicable to Purkinje fibre action and pacemaker potentials,” The Journal of Physiology, vol. 160, no. 2, pp. 317–352, 1962.

K. H. W. J. Ten Tusscher, D. Noble, P. J. Noble, and A. V. Panfilov, “A model for human ventricular tissue,” American Journal of Physiology-Heart and Circulatory Physiology, vol. 286, no. 4, pp. H1573–H1589, 2004.

C. H. Luo and Y. Rudy, “A dynamic model of the cardiac ventricular action potential. I. Simulations of ionic currents and concentration changes,” Circulation Research, vol. 74, no. 6, pp. 1071–1096, 1994.

M. Courtemanche, R. J. Ramirez, and S. Nattel, “Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model,” American Journal of Physiology-Heart and Circulatory Physiology, vol. 275, no. 1, pp. H301–H321, 1998.

R. FitzHugh, “Impulses and physiological states in theoretical models of nerve membrane,” Biophysical Journal, vol. 1, no. 6, pp. 445–466, 1961.

J. Nagumo, S. Arimoto, and S. Yoshizawa, “An active pulse transmission line simulating nerve axon,” Proceedings of the IRE, vol. 50, no. 10, pp. 2061–2070, 1962.

F. Fenton and A. Karma, “Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: filament instability and fibrillation,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 8, no. 1, pp. 20–47, 1998.

R. A. Oliver and W. Krassowska, “Reproducing cardiac restitution properties using the fenton-karma membrane model,” Annals of Biomedical Engineering, vol. 33, no. 7, pp. 907–911, 2005.

A. Bueno-Orovio, E. M. Cherry, and F. H. Fenton, “Minimal model for human ventricular action potentials in tissue,” Journal of Theoretical Biology, vol. 253, no. 3, pp. 544–560, 2008.

F. H. Fenton, S. J. Evans, and H. M. Hastings, “Memory in an excitable medium: a mechanism for spiral wave breakup in the low-excitability limit,” Physical Review Letters, vol. 83, no. 19, pp. 3964–3967, 1999.

S. E. Geneser, R. M. Kirby, and F. B. Sachse, “Sensitivity analysis of cardiac electrophysiological models using polynomial chaos,” in Proceedings of the 2005 IEEE Engineering in Medicine and Biology 27th Annual Conference, pp. 4042–4045, IEEE, New York, NY, USA, 2006.

Z. Syed, E. Vigmond, S. Nattel, and L. J. Leon, “Atrial cell action potential parameter fitting using genetic algorithms,” Medical & Biological Engineering & Computing, vol. 43, no. 5, pp. 561–571, 2005.

C. T. Bot, A. R. Kherlopian, F. A. Ortega, D. J. Christini, and T. Krogh-Madsen, “Rapid genetic algorithm optimization of a mouse computational model: benefits for anthropomorphization of neonatal mouse cardiomyocytes,” Frontiers in Physiology, vol. 3, p. 421, 2012.

W. Groenendaal, F. A. Ortega, A. R. Kherlopian, A. C. Zygmunt, T. Krogh-Madsen, and D. J. Christini, “Cell-specific cardiac electrophysiology models,” PLoS Computational Biology, vol. 11, no. 4, Article ID e1004242, 2015.

G. Duckett and D. Barkley, “Modeling the dynamics of cardiac action potentials,” Physical Review Letters, vol. 85, no. 4, pp. 884–887, 2000.

E. G. Tolkacheva, D. G. Schaeffer, D. J. Gauthier, and C. C. Mitchell, “Analysis of the Fenton-Karma model through an approximation by a one-dimensional map,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 12, no. 4, pp. 1034–1042, 2002.

T. Guo, A. Al Abed, N. H. Lovell, and S. Dokos, “Optimisation of a generic ionic model of cardiac myocyte electrical activity,” Computational Mathematical Methods in Medicine, vol. 2013, Article ID 706195, 20 pages, 2013.

S. Dokos and N. H. Lovell, “Parameter estimation in cardiac ionic models,” Progress in Biophysics Molecular Biology, vol. 85, no. 2-3, pp. 407–431, 2004.

A. X. Sarkar and E. A. Sobie, “Regression analysis for constraining free parameters in electrophysiological models of cardiac cells,” PLoS Computational Biology, vol. 6, no. 9, Article ID e1000914, 2010.

E. Angelino and M. P. Brenner, “Excitability constraints on voltage-gated sodium channels,” PLoS Computational Biology, vol. 3, no. 9, p. e177, 2007.

V. Biktashev and R. Suckley, “Non-Tikhonov asymptotic properties of cardiac excitability,” Physical Review Letters, vol. 93, no. 16, Article ID 168103, 2004.

D. Barkley, “A model for fast computer simulation of waves in excitable media,” Physica D: Nonlinear Phenomena, vol. 49, no. 1-2, pp. 61–70, 1991.

E. M. Izhikevich, Dynamical Systems in Neuroscience, MIT press, New York, NY, USA, 2007.

R. D. Simitev and V. N. Biktashev, “Asymptotics of conduction velocity restitution in models of electrical excitation in the heart,” Bulletin of Mathematical Biology, vol. 73, no. 1, pp. 72–115, 2011.

C. Mitchell and D. G. Schaeffer, “A two-current model for the dynamics of cardiac membrane,” Bulletin of Mathematical Biology, vol. 65, no. 5, pp. 767–793, 2003.

S. H. Sabzpoushan and A. Ghajarjazy, “A method for reduction of human ventricular action potential model,” Mathematical and Computer Modelling of Dynamical Systems, vol. 26, no. 1, pp. 1–30, 2020.

J. Marín-García, “Ion signaling and electrophysiological function,” in Signaling in the Heart, pp. 87–97, Springer, Berlin, Germany, 2011.

F. Burton and S. M. Cobbe, “Dispersion of ventricular repolarization and refractory period,” Cardiovascular Research, vol. 50, no. 1, pp. 10–23, 2001.

A. Garny and P. J. Hunter, “OpenCOR: a modular and interoperable approach to computational biology,” Frontiers in Physiology, vol. 6, p. 26, 2015.

C. M. Lloyd, J. R. Lawson, P. J. Hunter, and P. F. Nielsen, “The CellML model repository,” Bioinformatics, vol. 24, no. 18, pp. 2122-2123, 2008.

R. Noble, “Remodelling of calcium dynamics in guinea-pig ventricular cells,” 2001.

R. Courtemanche and J. Nattel, “Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model,” 1998.