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

Authors: S. H. Sabzpoushan; A. Ghajarjazy
Journal: Mathematical Problems in Engineering
 (2020)

---

## 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.


---

## Body

## 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].Figure1(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.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).
(a)(b)(c)(d)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, v, 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 vr, vt1,vp, vt2, T, and TE. Here, vr stands for the rest state or rest voltage, vt1 is the upstroke threshold, i.e., if the applied stimulus to the EXC puts the membrane voltage below vt1, the EXC will not be excited, and the membrane voltage, v, will return to vr, 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, vp; and then, after a time evolution TE (with the duration of T), it returns to the rest state, vr, through vt2. Here, vt2 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 Figure1, 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, vr. Phase 0 or depolarization phase is the counterpart of the fast upstroke phase in Figure 1(b), in which excitability threshold voltage (vt1) should be determined experimentally. In phase 1, a spike is generated in the VAP so that the spike voltage can be treated as vp 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), vt2 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 vp 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 T 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 Figure1(d) and general formalism of HH, our observation from a VTC is a voltage which is described by the following differential equation:(1)v˙=−1CIion1+Iion2+⋯+Iionk.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 notationIv for this group and 2-those which may be treated as functions of voltage and time simultaneously, with the symbolic notation Iv,t. In our simple model, we use one delegate from each group with symbols Ifv and Isv,t, respectively. Also, we refer to them as If and Is for simplicity. So, we can write(2)v˙=−1CIf+Is,where C stands for the membrane capacitance and IfandIs represent two types of ionic currents, sometimes named as fast and slow, respectively. The explicit expressions for these currents are given by(3)If=gfvv−Ef,(4)Is=gsv,tv−Es.In equation (3), gfv is a voltage-dependent membrane conductance, and Ef denotes a constant voltage (Nernst voltage by some ways). In equation (4), gsv,t is a voltage- and time-dependent (dynamic) membrane conductance, and Es 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(5)gsv,t=gs⋅nv,t,where gs is a constant conductance and nv,t delegates a voltage-dependent dynamic gating variable. The following equation describes the nonlinear differential equation of this gating variable:(6)n˙v,t=1τnv,∞−nv,t,where nv,∞, the steady state value of n at voltage v, is considered as a sigmoid function of v as follows:(7)nv,∞=11+expv12n−V/kn.In order to prepare more adaptability of our model with electrophysiology of VTCs, we decomposegfv as gf⋅mv in which gf is a constant conductance, and mv is a voltage-dependent gating variable. We also consider a sigmoid function of v for mv as(8)mv=11+expv12m−v/km.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, C,gf,v12m,km,Ef,gs,Es,τn,v12n,andkn, should be estimated so that they can simulate a real AP of a single VTC. It is desired that the last three parameters, τn,v12n,andkn, which are related to TE or T 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, vr,vt1,vp,vt2,andT, are fulfilled, and the effect of each parameter on these properties can be analysed explicitly.To present our novel graphical approach, we plotIfv and −Isv,t versus v on the same v-axis in Figure 2. Regarding equation (2), their intersection points mean the current balance or equilibrium points.Figure 2
A typical diagram of currentsIfv and −Isv,t versus v. On −∞<t≤0, two graphs have intersections at three points: A, B, and D. At t>0, the line of −Isv,t rotates clockwise around pivot point A. Therefore, B and D move toward each other and meet at vt2, where the line is tangent to the Ifv graph.CurrentIfv can be drawn easily as a function of v, i.e., versus v. According to equation (3), it intersects the v-axis at two points: v=Ef and where mv=0. Interestingly, the current −Isv,t at each time (tx) can be supposed as a line with the slope −gs⋅nv,tx that interceptsv-axis at v=Es. As we mentioned above, the possible intersection points of the two graphs at each tx indicate the equilibrium points of the system, i.e., v˙=0.Regarding the above discussion, as time increases, ifnv,t be about increasing, the line −Isv,t rotates clockwise around its pivot point, v=Es, and if nv,t be decreasing, the line −Isv,t rotates counterclockwise. We see that the pivot point v=Es (that we will set vr≈Es in future) is a key point.Now, suppose we have adjusted the model parameters att=0 so that Ifv and −Isv,0 have three intersection points (A, B, and D), where vA=vr, vB=vt1, and vD=vp as depicted in Figure 2. All of these points are equilibrium points because the equality condition is satisfied (Ifv=−Isv,0) at all of them. Assume that the model is at the rest state, i.e., v=vr; we apply a subthreshold excitation that increases v just a little (v=v0) so that v0<vt1 in this situation. It is clear from Figure 2 that Isv0,0>Ifv0; consequently, regarding equation (2), the membrane potential will be pushed back to the rest state, vr. 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, thenv increases to v1 so that v1>vt1. Regarding Figure 2 and equation (2), Isv1,0<Ifv1 leads to the increase in voltage. As a consequence, the membrane potential will be attracted by the next equilibrium point vp. All of these mean that the model will be excited by a superthreshold excitation.The above paragraphs reveal that the graphical illustration in Figure2 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 fromvr to vp, through vt1, and develops from vt2 to vr. A visual deduction from Figure 2 is(9)Es<vr<vt1<vt2<vp<Ef.Since time has not an important role in the excitability section, we can consider that all things happen instantaneously. Therefore,nv,t has not been activated throughout this section and remains at its value before the stimulated state, nvr,∞. Assuming C=1 for simplicity, we rewrite equations (2) and (6), respectively, as(10)v˙=−gf⋅mvv−Ef+gs⋅nvr,∞v−Es,(11)n˙v,t=0.CurrentIfv in equation (10) has two important roles. First, it makes a positive instantaneous feedback and a fast regenerative process for pushing v from vt1 to vp. Second, as we will see later, it opposes Isv,t in the plateau phase (v evolution from vp toward vt2) and controls the decrements of v as well. Concerning the excitability definition, the amplitude of Ifv should be small on the vr to vt1 interval to suppress the subthreshold stimulus and large on the vt1 to vp interval in order to cause fast regenerative process. The above required properties for If are provided by the definition of mv in equation (8).AssumingIfvr≈0, we can write equation (12) in the triangle with corners vr,vp, and D:(12)Ifvt1Ifvp=vt1−vrvp−vr.We define parameterα as the ratio of m-gate in two distinct voltages (vt1 and vp) as follows:(13)mvt1mvp=α,0<α<1.Substituting equation (13) into equation (3), one gets(14)gf⋅mvt1vt1−Ef=α⋅gf⋅mvpvt1−Ef.Using equation (14), we can write(15)Ifvt1Ifvp=αvt1−Efvp−Ef.Solving algebraic equations (12) and (15) for Ef, one gets(16)Ef=αvt1vp−vr−vpvt1−vrαvp−vr−vt1−vr.(16) means that, by having desired characteristics, vr,vt1, and vp, for a cell and choosing 0<α<1, we can get one of the model parameters, Ef.Considering−Isv,nvr,∞ line, at the vp fixed point, we have(17)Isvp,nvr,∞+Ifvp=0.For simplicity and without loss of generality, we can scale currents so that(18)Ifvp=−1.By substituting equations (4) and (18) into equation (17), we get(19)gsnvr,∞vp−Es−1=0.Therefore,(20)gsnvr,∞=1vp−Es.Writing current balance at another equilibrium point (vt1) leads to(21)Ifvt1=−Isvt1.Using equations (15), (18), and (20) in equation (21), we get(22)αIfvpvt1−Efvp−Ef=−gsnvr,∞vt1−Es.Solving equation (22) for α gives(23)α=vp−Efvt1−EfIsvt1=vp−Efvt1−Esvt1−Efvp−Es.After some algebraic manipulations on equation (23), we get equation (24), which is a counterpart of equation (16):(24)Ef=αvt1vp−Es−vpvt1−Esαvp−Es−vt1−Es.Satisfaction of equations (16) and (24) simultaneously results in the following equation:(25)Es=vr.It means that another parameter of the model,Es, is determined.Writing current balance at stable fixed pointvr leads to Ifvr=Isvr,nvr,∞=0; therefore, we can treat α as a free parameter and choose it small enough to satisfy Ifvr≈0.According to them-gate definition, if mvt1<1, then mvr<mvt1 which leads to mvr≪1, so the hope is that the rough condition Ifvr≈0 be satisfied better. In this condition,(26)mvt1=11+expv12m−vt1/km=α,α≪1.By some algebraic manipulations on equation (26), we get a relation between parameters of m-gate:(27)km=1βv12m−vt1,β=ln1−αα.In addition, the instantaneous process on thevt2 to vr interval can be convincing by defining vt2 as a saddle fixed point at time T. It means that −Isv,T should be tangent to Ifv at point vt2, which is fulfilled by the following condition:(28)1km.gfmvt21−mvt2vt2−Ef+gfmvt2=Ifvt2vt2−vr.By substituting equation (27) into equation (28), we get(29)11+exp−βv12m−vt2/v12m−vt1=1β.vr−Efvt2−Ef−vt2−vrv12m−vt1.Solving equation (29) leads to the determination of v12m, 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 v12m, and the cross section point is the solution. Afterward, using obtained v12m in equation (27) results in the determination of km. Now, having v12m and km, values of m-gate can be found at any voltage. Therefore, gf can be determined by the following equation:(30)gf=−1+expv12m−vp/kmvp−Ef.In conclusion, we see that the parametersv12m, km, gf, Es, and Ef 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 ofIs about its pivot point v=vr. This rotation is attributed to gradual activation of nv,t, which causes a time course of voltage progress from vp to vt2. It is noteworthy that although nv,t 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 assumenv,t as an activation gate and choose an arbitrary small value for it at the rest state:(31)nvr,∞=11+expv12n−v/kn=γ,γ≪1.Using equation (31), satisfying gsnvr,∞ in equation (20), we can determine gs:(32)gs=gsnvr,∞γ.Recalling our graphical approach, the line−Isv,t should be tangent to the Ifv graph at v=vt2,t=T, so(33)−Isvt2,T=Ifvt2.Then, by considering equations (3) and (4), one gets(34)nvt2,T=−Ifvt2gsvt2−Es=θ.Having two values ofn, in equations (31) and (34), we can determine v12n and kn as follows:(35)kn=vt2−vrLn1/γ−1/1/θ−1,(36)v12n=vt2+knLn1/θ−1.Having in mind the above discussions about excitability and TE sections and considering Figure2, as well as presented equations, we implemented numeric calculations of our graphical approach in Matlab as follows.At the first time stepΔt after excitation, nvr,∞raises toward nvp,∞; this increment makes a clockwise rotation of the −Isv,t line, i.e., generates −Isv+Δv,t+Δt and forms a new fixed point (e.g., point d shown in Figure 2). This new fixed point creates new n˙, as well. As time progresses, a sequence of n˙v,t and consequently nv,t are generated. To satisfy TE of wanted AP and inducing fast repolarization at vt2, the sequence of the generated two fixed points (d, b) should meet each other at time t=T and at voltage vt2. When this happens, the regenerative process occurs to get the voltage back to vr quickly. Regarding equation (6), the time step Δt should be small enough so that we can capture the dynamics of nv,t during the rotation of the −Isv,t line. In Figure 3, the general algorithm for the determination of the successive intersection points of the −Isv,t line with the Ifv graph is summarized.Figure 3
The algorithm for the determination of successive intersection points of the−Isv,t line with the Ifv graph.In Figure4, 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.


## 2.1. Graphical Approach and Model Equations

Regarding Figure1(d) and general formalism of HH, our observation from a VTC is a voltage which is described by the following differential equation:(1)v˙=−1CIion1+Iion2+⋯+Iionk.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 notationIv for this group and 2-those which may be treated as functions of voltage and time simultaneously, with the symbolic notation Iv,t. In our simple model, we use one delegate from each group with symbols Ifv and Isv,t, respectively. Also, we refer to them as If and Is for simplicity. So, we can write(2)v˙=−1CIf+Is,where C stands for the membrane capacitance and IfandIs represent two types of ionic currents, sometimes named as fast and slow, respectively. The explicit expressions for these currents are given by(3)If=gfvv−Ef,(4)Is=gsv,tv−Es.In equation (3), gfv is a voltage-dependent membrane conductance, and Ef denotes a constant voltage (Nernst voltage by some ways). In equation (4), gsv,t is a voltage- and time-dependent (dynamic) membrane conductance, and Es 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(5)gsv,t=gs⋅nv,t,where gs is a constant conductance and nv,t delegates a voltage-dependent dynamic gating variable. The following equation describes the nonlinear differential equation of this gating variable:(6)n˙v,t=1τnv,∞−nv,t,where nv,∞, the steady state value of n at voltage v, is considered as a sigmoid function of v as follows:(7)nv,∞=11+expv12n−V/kn.In order to prepare more adaptability of our model with electrophysiology of VTCs, we decomposegfv as gf⋅mv in which gf is a constant conductance, and mv is a voltage-dependent gating variable. We also consider a sigmoid function of v for mv as(8)mv=11+expv12m−v/km.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, C,gf,v12m,km,Ef,gs,Es,τn,v12n,andkn, should be estimated so that they can simulate a real AP of a single VTC. It is desired that the last three parameters, τn,v12n,andkn, which are related to TE or T 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, vr,vt1,vp,vt2,andT, are fulfilled, and the effect of each parameter on these properties can be analysed explicitly.To present our novel graphical approach, we plotIfv and −Isv,t versus v on the same v-axis in Figure 2. Regarding equation (2), their intersection points mean the current balance or equilibrium points.Figure 2
A typical diagram of currentsIfv and −Isv,t versus v. On −∞<t≤0, two graphs have intersections at three points: A, B, and D. At t>0, the line of −Isv,t rotates clockwise around pivot point A. Therefore, B and D move toward each other and meet at vt2, where the line is tangent to the Ifv graph.CurrentIfv can be drawn easily as a function of v, i.e., versus v. According to equation (3), it intersects the v-axis at two points: v=Ef and where mv=0. Interestingly, the current −Isv,t at each time (tx) can be supposed as a line with the slope −gs⋅nv,tx that interceptsv-axis at v=Es. As we mentioned above, the possible intersection points of the two graphs at each tx indicate the equilibrium points of the system, i.e., v˙=0.Regarding the above discussion, as time increases, ifnv,t be about increasing, the line −Isv,t rotates clockwise around its pivot point, v=Es, and if nv,t be decreasing, the line −Isv,t rotates counterclockwise. We see that the pivot point v=Es (that we will set vr≈Es in future) is a key point.Now, suppose we have adjusted the model parameters att=0 so that Ifv and −Isv,0 have three intersection points (A, B, and D), where vA=vr, vB=vt1, and vD=vp as depicted in Figure 2. All of these points are equilibrium points because the equality condition is satisfied (Ifv=−Isv,0) at all of them. Assume that the model is at the rest state, i.e., v=vr; we apply a subthreshold excitation that increases v just a little (v=v0) so that v0<vt1 in this situation. It is clear from Figure 2 that Isv0,0>Ifv0; consequently, regarding equation (2), the membrane potential will be pushed back to the rest state, vr. 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, thenv increases to v1 so that v1>vt1. Regarding Figure 2 and equation (2), Isv1,0<Ifv1 leads to the increase in voltage. As a consequence, the membrane potential will be attracted by the next equilibrium point vp. All of these mean that the model will be excited by a superthreshold excitation.The above paragraphs reveal that the graphical illustration in Figure2 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 fromvr to vp, through vt1, and develops from vt2 to vr. A visual deduction from Figure 2 is(9)Es<vr<vt1<vt2<vp<Ef.Since time has not an important role in the excitability section, we can consider that all things happen instantaneously. Therefore,nv,t has not been activated throughout this section and remains at its value before the stimulated state, nvr,∞. Assuming C=1 for simplicity, we rewrite equations (2) and (6), respectively, as(10)v˙=−gf⋅mvv−Ef+gs⋅nvr,∞v−Es,(11)n˙v,t=0.CurrentIfv in equation (10) has two important roles. First, it makes a positive instantaneous feedback and a fast regenerative process for pushing v from vt1 to vp. Second, as we will see later, it opposes Isv,t in the plateau phase (v evolution from vp toward vt2) and controls the decrements of v as well. Concerning the excitability definition, the amplitude of Ifv should be small on the vr to vt1 interval to suppress the subthreshold stimulus and large on the vt1 to vp interval in order to cause fast regenerative process. The above required properties for If are provided by the definition of mv in equation (8).AssumingIfvr≈0, we can write equation (12) in the triangle with corners vr,vp, and D:(12)Ifvt1Ifvp=vt1−vrvp−vr.We define parameterα as the ratio of m-gate in two distinct voltages (vt1 and vp) as follows:(13)mvt1mvp=α,0<α<1.Substituting equation (13) into equation (3), one gets(14)gf⋅mvt1vt1−Ef=α⋅gf⋅mvpvt1−Ef.Using equation (14), we can write(15)Ifvt1Ifvp=αvt1−Efvp−Ef.Solving algebraic equations (12) and (15) for Ef, one gets(16)Ef=αvt1vp−vr−vpvt1−vrαvp−vr−vt1−vr.(16) means that, by having desired characteristics, vr,vt1, and vp, for a cell and choosing 0<α<1, we can get one of the model parameters, Ef.Considering−Isv,nvr,∞ line, at the vp fixed point, we have(17)Isvp,nvr,∞+Ifvp=0.For simplicity and without loss of generality, we can scale currents so that(18)Ifvp=−1.By substituting equations (4) and (18) into equation (17), we get(19)gsnvr,∞vp−Es−1=0.Therefore,(20)gsnvr,∞=1vp−Es.Writing current balance at another equilibrium point (vt1) leads to(21)Ifvt1=−Isvt1.Using equations (15), (18), and (20) in equation (21), we get(22)αIfvpvt1−Efvp−Ef=−gsnvr,∞vt1−Es.Solving equation (22) for α gives(23)α=vp−Efvt1−EfIsvt1=vp−Efvt1−Esvt1−Efvp−Es.After some algebraic manipulations on equation (23), we get equation (24), which is a counterpart of equation (16):(24)Ef=αvt1vp−Es−vpvt1−Esαvp−Es−vt1−Es.Satisfaction of equations (16) and (24) simultaneously results in the following equation:(25)Es=vr.It means that another parameter of the model,Es, is determined.Writing current balance at stable fixed pointvr leads to Ifvr=Isvr,nvr,∞=0; therefore, we can treat α as a free parameter and choose it small enough to satisfy Ifvr≈0.According to them-gate definition, if mvt1<1, then mvr<mvt1 which leads to mvr≪1, so the hope is that the rough condition Ifvr≈0 be satisfied better. In this condition,(26)mvt1=11+expv12m−vt1/km=α,α≪1.By some algebraic manipulations on equation (26), we get a relation between parameters of m-gate:(27)km=1βv12m−vt1,β=ln1−αα.In addition, the instantaneous process on thevt2 to vr interval can be convincing by defining vt2 as a saddle fixed point at time T. It means that −Isv,T should be tangent to Ifv at point vt2, which is fulfilled by the following condition:(28)1km.gfmvt21−mvt2vt2−Ef+gfmvt2=Ifvt2vt2−vr.By substituting equation (27) into equation (28), we get(29)11+exp−βv12m−vt2/v12m−vt1=1β.vr−Efvt2−Ef−vt2−vrv12m−vt1.Solving equation (29) leads to the determination of v12m, 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 v12m, and the cross section point is the solution. Afterward, using obtained v12m in equation (27) results in the determination of km. Now, having v12m and km, values of m-gate can be found at any voltage. Therefore, gf can be determined by the following equation:(30)gf=−1+expv12m−vp/kmvp−Ef.In conclusion, we see that the parametersv12m, km, gf, Es, and Ef 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 ofIs about its pivot point v=vr. This rotation is attributed to gradual activation of nv,t, which causes a time course of voltage progress from vp to vt2. It is noteworthy that although nv,t 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 assumenv,t as an activation gate and choose an arbitrary small value for it at the rest state:(31)nvr,∞=11+expv12n−v/kn=γ,γ≪1.Using equation (31), satisfying gsnvr,∞ in equation (20), we can determine gs:(32)gs=gsnvr,∞γ.Recalling our graphical approach, the line−Isv,t should be tangent to the Ifv graph at v=vt2,t=T, so(33)−Isvt2,T=Ifvt2.Then, by considering equations (3) and (4), one gets(34)nvt2,T=−Ifvt2gsvt2−Es=θ.Having two values ofn, in equations (31) and (34), we can determine v12n and kn as follows:(35)kn=vt2−vrLn1/γ−1/1/θ−1,(36)v12n=vt2+knLn1/θ−1.Having in mind the above discussions about excitability and TE sections and considering Figure2, as well as presented equations, we implemented numeric calculations of our graphical approach in Matlab as follows.At the first time stepΔt after excitation, nvr,∞raises toward nvp,∞; this increment makes a clockwise rotation of the −Isv,t line, i.e., generates −Isv+Δv,t+Δt and forms a new fixed point (e.g., point d shown in Figure 2). This new fixed point creates new n˙, as well. As time progresses, a sequence of n˙v,t and consequently nv,t are generated. To satisfy TE of wanted AP and inducing fast repolarization at vt2, the sequence of the generated two fixed points (d, b) should meet each other at time t=T and at voltage vt2. When this happens, the regenerative process occurs to get the voltage back to vr quickly. Regarding equation (6), the time step Δt should be small enough so that we can capture the dynamics of nv,t during the rotation of the −Isv,t line. In Figure 3, the general algorithm for the determination of the successive intersection points of the −Isv,t line with the Ifv graph is summarized.Figure 3
The algorithm for the determination of successive intersection points of the−Isv,t line with the Ifv graph.In Figure4, 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.


## 2.1.1. Graphical Approach

In a simple model which is described by equations (2)–(8), ten parameters, C,gf,v12m,km,Ef,gs,Es,τn,v12n,andkn, should be estimated so that they can simulate a real AP of a single VTC. It is desired that the last three parameters, τn,v12n,andkn, which are related to TE or T 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, vr,vt1,vp,vt2,andT, are fulfilled, and the effect of each parameter on these properties can be analysed explicitly.To present our novel graphical approach, we plotIfv and −Isv,t versus v on the same v-axis in Figure 2. Regarding equation (2), their intersection points mean the current balance or equilibrium points.Figure 2
A typical diagram of currentsIfv and −Isv,t versus v. On −∞<t≤0, two graphs have intersections at three points: A, B, and D. At t>0, the line of −Isv,t rotates clockwise around pivot point A. Therefore, B and D move toward each other and meet at vt2, where the line is tangent to the Ifv graph.CurrentIfv can be drawn easily as a function of v, i.e., versus v. According to equation (3), it intersects the v-axis at two points: v=Ef and where mv=0. Interestingly, the current −Isv,t at each time (tx) can be supposed as a line with the slope −gs⋅nv,tx that interceptsv-axis at v=Es. As we mentioned above, the possible intersection points of the two graphs at each tx indicate the equilibrium points of the system, i.e., v˙=0.Regarding the above discussion, as time increases, ifnv,t be about increasing, the line −Isv,t rotates clockwise around its pivot point, v=Es, and if nv,t be decreasing, the line −Isv,t rotates counterclockwise. We see that the pivot point v=Es (that we will set vr≈Es in future) is a key point.Now, suppose we have adjusted the model parameters att=0 so that Ifv and −Isv,0 have three intersection points (A, B, and D), where vA=vr, vB=vt1, and vD=vp as depicted in Figure 2. All of these points are equilibrium points because the equality condition is satisfied (Ifv=−Isv,0) at all of them. Assume that the model is at the rest state, i.e., v=vr; we apply a subthreshold excitation that increases v just a little (v=v0) so that v0<vt1 in this situation. It is clear from Figure 2 that Isv0,0>Ifv0; consequently, regarding equation (2), the membrane potential will be pushed back to the rest state, vr. 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, thenv increases to v1 so that v1>vt1. Regarding Figure 2 and equation (2), Isv1,0<Ifv1 leads to the increase in voltage. As a consequence, the membrane potential will be attracted by the next equilibrium point vp. All of these mean that the model will be excited by a superthreshold excitation.The above paragraphs reveal that the graphical illustration in Figure2 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 fromvr to vp, through vt1, and develops from vt2 to vr. A visual deduction from Figure 2 is(9)Es<vr<vt1<vt2<vp<Ef.Since time has not an important role in the excitability section, we can consider that all things happen instantaneously. Therefore,nv,t has not been activated throughout this section and remains at its value before the stimulated state, nvr,∞. Assuming C=1 for simplicity, we rewrite equations (2) and (6), respectively, as(10)v˙=−gf⋅mvv−Ef+gs⋅nvr,∞v−Es,(11)n˙v,t=0.CurrentIfv in equation (10) has two important roles. First, it makes a positive instantaneous feedback and a fast regenerative process for pushing v from vt1 to vp. Second, as we will see later, it opposes Isv,t in the plateau phase (v evolution from vp toward vt2) and controls the decrements of v as well. Concerning the excitability definition, the amplitude of Ifv should be small on the vr to vt1 interval to suppress the subthreshold stimulus and large on the vt1 to vp interval in order to cause fast regenerative process. The above required properties for If are provided by the definition of mv in equation (8).AssumingIfvr≈0, we can write equation (12) in the triangle with corners vr,vp, and D:(12)Ifvt1Ifvp=vt1−vrvp−vr.We define parameterα as the ratio of m-gate in two distinct voltages (vt1 and vp) as follows:(13)mvt1mvp=α,0<α<1.Substituting equation (13) into equation (3), one gets(14)gf⋅mvt1vt1−Ef=α⋅gf⋅mvpvt1−Ef.Using equation (14), we can write(15)Ifvt1Ifvp=αvt1−Efvp−Ef.Solving algebraic equations (12) and (15) for Ef, one gets(16)Ef=αvt1vp−vr−vpvt1−vrαvp−vr−vt1−vr.(16) means that, by having desired characteristics, vr,vt1, and vp, for a cell and choosing 0<α<1, we can get one of the model parameters, Ef.Considering−Isv,nvr,∞ line, at the vp fixed point, we have(17)Isvp,nvr,∞+Ifvp=0.For simplicity and without loss of generality, we can scale currents so that(18)Ifvp=−1.By substituting equations (4) and (18) into equation (17), we get(19)gsnvr,∞vp−Es−1=0.Therefore,(20)gsnvr,∞=1vp−Es.Writing current balance at another equilibrium point (vt1) leads to(21)Ifvt1=−Isvt1.Using equations (15), (18), and (20) in equation (21), we get(22)αIfvpvt1−Efvp−Ef=−gsnvr,∞vt1−Es.Solving equation (22) for α gives(23)α=vp−Efvt1−EfIsvt1=vp−Efvt1−Esvt1−Efvp−Es.After some algebraic manipulations on equation (23), we get equation (24), which is a counterpart of equation (16):(24)Ef=αvt1vp−Es−vpvt1−Esαvp−Es−vt1−Es.Satisfaction of equations (16) and (24) simultaneously results in the following equation:(25)Es=vr.It means that another parameter of the model,Es, is determined.Writing current balance at stable fixed pointvr leads to Ifvr=Isvr,nvr,∞=0; therefore, we can treat α as a free parameter and choose it small enough to satisfy Ifvr≈0.According to them-gate definition, if mvt1<1, then mvr<mvt1 which leads to mvr≪1, so the hope is that the rough condition Ifvr≈0 be satisfied better. In this condition,(26)mvt1=11+expv12m−vt1/km=α,α≪1.By some algebraic manipulations on equation (26), we get a relation between parameters of m-gate:(27)km=1βv12m−vt1,β=ln1−αα.In addition, the instantaneous process on thevt2 to vr interval can be convincing by defining vt2 as a saddle fixed point at time T. It means that −Isv,T should be tangent to Ifv at point vt2, which is fulfilled by the following condition:(28)1km.gfmvt21−mvt2vt2−Ef+gfmvt2=Ifvt2vt2−vr.By substituting equation (27) into equation (28), we get(29)11+exp−βv12m−vt2/v12m−vt1=1β.vr−Efvt2−Ef−vt2−vrv12m−vt1.Solving equation (29) leads to the determination of v12m, 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 v12m, and the cross section point is the solution. Afterward, using obtained v12m in equation (27) results in the determination of km. Now, having v12m and km, values of m-gate can be found at any voltage. Therefore, gf can be determined by the following equation:(30)gf=−1+expv12m−vp/kmvp−Ef.In conclusion, we see that the parametersv12m, km, gf, Es, and Ef 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 ofIs about its pivot point v=vr. This rotation is attributed to gradual activation of nv,t, which causes a time course of voltage progress from vp to vt2. It is noteworthy that although nv,t 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 assumenv,t as an activation gate and choose an arbitrary small value for it at the rest state:(31)nvr,∞=11+expv12n−v/kn=γ,γ≪1.Using equation (31), satisfying gsnvr,∞ in equation (20), we can determine gs:(32)gs=gsnvr,∞γ.Recalling our graphical approach, the line−Isv,t should be tangent to the Ifv graph at v=vt2,t=T, so(33)−Isvt2,T=Ifvt2.Then, by considering equations (3) and (4), one gets(34)nvt2,T=−Ifvt2gsvt2−Es=θ.Having two values ofn, in equations (31) and (34), we can determine v12n and kn as follows:(35)kn=vt2−vrLn1/γ−1/1/θ−1,(36)v12n=vt2+knLn1/θ−1.Having in mind the above discussions about excitability and TE sections and considering Figure2, as well as presented equations, we implemented numeric calculations of our graphical approach in Matlab as follows.At the first time stepΔt after excitation, nvr,∞raises toward nvp,∞; this increment makes a clockwise rotation of the −Isv,t line, i.e., generates −Isv+Δv,t+Δt and forms a new fixed point (e.g., point d shown in Figure 2). This new fixed point creates new n˙, as well. As time progresses, a sequence of n˙v,t and consequently nv,t are generated. To satisfy TE of wanted AP and inducing fast repolarization at vt2, the sequence of the generated two fixed points (d, b) should meet each other at time t=T and at voltage vt2. When this happens, the regenerative process occurs to get the voltage back to vr quickly. Regarding equation (6), the time step Δt should be small enough so that we can capture the dynamics of nv,t during the rotation of the −Isv,t line. In Figure 3, the general algorithm for the determination of the successive intersection points of the −Isv,t line with the Ifv graph is summarized.Figure 3
The algorithm for the determination of successive intersection points of the−Isv,t line with the Ifv graph.In Figure4, 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.Table1 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 −Isv,t line with the Ifv curve, and the cross sections at 50, −28, −49, and −92.8 mv voltages correspond to vp, vt2,vt1, and vr, 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.1gfigsongsoEfiEsokmv12mknv12nτγ=0.10.02090.0070.0797.89−92.84.96−38.128.8−29.4375.5η=0.001Figure 5
Reproducing AP excitability generated by the Noble-2001 model. (a) Proposed graphical approach. (b) Simulated AP.
(a)(b)Table 2
Comparing the excitability features of the AP simulated by the proposed model with the ones generated by the Noble-2001 model.
vrvt1vt2vpTτIfivrNoble-2001−92.8 mv−49 mv−28 mv50 mv160 msec——Parametersα=0.1−92.7−49−284916075.43−0.000064α=0.2−92.55−48.91−284916080.52−0.000349In 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.1gfigsongsoEfiEsokmv12mknv12nτγ=0.10.01600.00950.187.0659−81.194.1416−47.780021.3548−34.2686135η=0.001Figure 6
Reproducing excitability of the AP generated by the CRN-1998 model. (a) Proposed graphical approach. (b) Simulated AP.
(a)(b)Table 4
Comparing the excitability features of the AP simulated by the proposed model with the ones generated by the CRN model.
vrvt1vt2vpTCRN−81.19 mv−56.88 mv−40 mv24.49 mv215.6 msecParametersα=0.1−81.19−56−4023.8214.3Reproduced 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 Figure1 formalism.In our paper, the focus of attention was on excitability simulation, i.e., excitability section of Figures1(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 Figures5(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 fornv,t, but our approach is flexible so that we can consider different speed profiles for nv,t. To be more clear, we may even attribute higher-order dynamics to nv,t in Figure 2, i.e., rotation of −Isv,t 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 −Isv,t 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 n,n˙,n¨ 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 n as depicted in Figures 8(b)–8(d).Figure 7
Using higher-order dynamics forn 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) n time course profile. (c) n˙ speed profile. (d) n¨ acceleration profile.
(a)(b)(c)(d)Figure 8
Using sample second-order dynamics forn 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) n time course profile. (c) n˙ speed profile. (d) n¨ acceleration profile. The parameter values of the excitability section are the same as Table 1.
(a)(b)(c)(d)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 liken4, m3h,m3hj, nh, … in the electrophysiology literature [3, 5, 12, 26].In these cases, we can assume a term like nh 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 ofnv,t 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 −Isv,t line after leaving tangency at the vt2 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.



---

*Source: 3704523-2020-07-22.xml*