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)
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 writewhere 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 writewhere 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.
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.
In Figure 4, the whole steps and sequences of parameter estimation of our approach are illustrated.
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.
(a)
(b)
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.
(a)
(b)
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)
(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 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)