A two-variable model robust to pacemaker behaviour for the dynamics of the cardiac action potential

Highlights • Simplified cardiomyocyte electrophysiology model optimised for patient specific modelling.• Modified Mitchell and Schaeffer model.• No spurious pacemaker behaviour.• Suitable for parameter fitting.• Consistent leading order APD approximation.


Introduction
Cardiac ionic models are mathematical models describing the electrical response of a cardiac myocyte following an applied electrical stimulus.
When an electrical stimulus is applied to a cardiac myocyte, an action potential is generated by the flux of ionic species across the cell membrane. Complex mathematical models [2,15,24] describe the ionic current generated by each ionic species. Though physiologically accurate, these models are expensive to solve numerically due to their large number of state variables and their non-linear formulation; moreover, when personalising such models, an additional challenge arises from the large number of parameters to constrain.
In contrast, phenomenological ionic models aim to describe the collective effects of the ionic currents by a smaller number of state variables (usually 2-3) and parameters, and are obtained either by simplifying a complex ionic model [7,16] , or by trying to reproduce the shape of the action potential [1,8] . In the field of personalised models, particular emphasis has been given to the models with two state variables, due to their numerical tractability and the reduced number of parameters to be constrained [4,10,20,23,27] .
Among the available two state-variable ionic models, the one introduced by Mitchell and Schaeffer (MS) [16] is often employed in ventricular electro-physiology inverse problems [4,21,23] . This model is capable of reproducing the shape of the action potential and the restitution properties of the action potential duration (APD), [12,14] . The MS model is characterised by 5 parameters. It is obtained by simplifying the model proposed by Fenton and Karma, [7] (3 state variables, 13 parameters), which in turn is obtained by simplifying the more complex Luo Rudy I [15] biophysical ionic model (8 state variables, 63 parameters).
For some choices of the parameters, the MS model suffers from the so called "pacemaker cell behaviour". That is, the transmembrane potential cyclically depolarises and repolarises in the absence of any applied external stimulus, as depicted in Fig. 1 A. In Fig. 1 B (blue line) the solution of the MS model under pacemaker behaviour is depicted in the phase space. From an analysis of the MS model equations, it is possible to analytically define the relation between the gate variable and the trans-membrane potential delimiting the values of the state variables producing a depolarisation. This curve is called a nullcline, and it is split into a left branch, ( v − m ) and a right branch, ( v + m ) depicted in Fig. 1 B. Once activated when the system moves towards the initial (at rest) condition, if the phase portrait crosses the nullcline branch v − m , the system "falls" into a condition where a depolarisation occurs, and a new action potential will be produced even though no external stimuli were applied.  pacemaker cell behaviour do not delimit a closed region of the parameter space. This phenomenon has been observed and reported both in 0D and tissue models, [5,22] . However, to the best of our knowledge there are no criteria for determining a priori which combination of parameters will produce pacemaker activity.
For patient specific modelling, this unwanted behaviour represents a problem for parameter estimation, since it necessitates a stability test for each estimated parameter set, [5] .
Particularly affected by this phenomenon are the sequential data assimilation techniques [3,9,11] : since the values of the parameters are sequentially updated to minimise the discrepancy between the output of the model and the measurements of the system under study, if a combination of parameters yielding pacemaker behaviour is produced, the algorithm adopted could become unstable and diverge.
To overcome these limitations, in this paper we derive and introduce a two state-variable ionic model that describes the action potential phases by 5 parameters, has the same benefits as the MS model, and is robust to pacemaker behaviour. These characteristics make it suitable for generating personalised electrophysiology models for clinical applications, in particular when a sequential data assimilation technique is employed. This paper is organised as follows: in Section 2 we introduce the mathematical formulation of the new ionic model; in Section 3 we prove the absence of pacemaker behaviour in the phase plane; in Section 4 we derive an asymptotic derivation of the restitution curves and compare it to the one described in [16] for the standard MS model; and in Section 5 we compare the solutions and the restitution properties of the new model with the MS model for some numerical examples.

The modified Mitchell-Schaeffer ionic model
The standard MS ionic model [16] describes the ionic currents that flow across the cell membrane with a gated-inward ionic current, representing the current produced by the flux of the sodium ions, and an ungated outward ionic current, representing the current produced by the flux of the potassium ions. The two state variables characterising the MS model represent the electric potential of the cell membrane and the gate dynamics of the sodium ion channels. The model can be written in the generalised form proposed by [6] : where J stim is an externally applied electrical stimulus, v m is the trans-membrane potential, h is the gate variable of the inward current, v gate is the activation threshold potential and The effect of the potassium ion current is negligible when v m ≤ v gate and the cell is returning to a quiescent state. It is possible to introduce gating effects to the outward ionic current with the complement of the gate variable, (1 −h ), adopting an expression similar to the one introduced in [25] . This leads to the following system of ODEs: Remark. The robustness to pacemaker behaviour is obtained by modifying the cubic polynomial on the right hand side of Eq. (1) . The addition of gating on the outward current has three advantages: first, when h = 1 the threshold value of the transmembrane voltage above which an action potential is triggered is equal to v gate ; second, for h = 1 the transmembrane potential at the end of the upstroke will be v m = 1 ; third, as will be discussed in Section 4 , the analytical solution of the mMS model coincides to the analytical solution of the MS model for the same set of ionic parameters for a particular choice of v gate .

Robustness to pacemaker behaviour
The study of the robustness of the mMS model to pacemaker behaviour consists in determining if the condition ∂ v m / ∂ t > 0 could occur without an external stimulus being applied. In the phase plane, the nullclines of Eq. (1) for the MS model and Eq. (3) for the mMS model define the boundary between the region where the depolarisation occurs and the region where the system naturally evolves towards the rest condition; when h recovers towards its rest state h = 1 , if the phase portrait crosses the nullcline branch v − m , as depicted in Fig. 1 B, then the system depolarises again and a new action potential is produced without an external stimulus having been applied.
The robustness of the mMS model to pacemaker behaviour is shown by proving that its phase portrait cannot intersect the nullcline branch v − m unless an external stimulus is applied. To this aim, we first derive the nullcline of Eq.

Nullclines
Denoting the minimum value of h on the nullclines in the MS model by the nullclines of MS are defined in [16] as follows: denoting by the minimum value of h on the nullcline of the mMS model, the nullclines of the new mMS model described in Eq. (3) are defined as follows:  This point characterises the value of v m at the end of the depolarisation upstroke and it is equal to (1, 1) for mMS ( ) and to ( 1 . This point represents the minimum value an external stimulus has to rise v m , to produce an action potential when the system is fully recovered. It is equal to ( v gate , 1) ( ) for mMS and to In the phase plane, the phase portrait is defined as the curve described by the values ( v m ( t ), h ( t )) which constitute the solution of the ionic model. The phase portraits of MS ( ) and mMS ( ) are depicted in the right panel of Fig. 2 for a non-pacemaker solution; the phase portrait for a pacemaker solution for the MS model is depicted in Fig. 2 B.

Robustness
The nullcline and characteristic points were defined above. From these points we show that the phase portrait of mMS cannot cross the nullcline branch v − m during the recovery of h unless an external stimulus is applied. To simplify the notation, we denote the value of the trans-membrane potential of Point 2 in mMS by We then split the phase portrait into two parts: the repolarisation, defined by the interval v gate ≤ v m ≤ v * m for mMS and the interval v gate ≤ v m ≤ 1/2 for MS, and the recovery, when h returns to its initial value, defined by v m < v gate .
For both models, during the repolarisation v m approaches v gate at a rate proportional to 1/ τ out . Since v m ≥ v gate , h decreases at a rate proportional to 1/ τ close and the phase portrait moves away from the nullcline branch v − m , lying in the region of the phase plane where no depolarisation can occur.
When v m < v gate , h recovers at a rate proportional to 1/ τ open , while the time derivative of v m is still negative, as depicted in We tested the robustness of the MS and the mMS models to pacemaker behaviour on a regular grid set of 57,800 parameters chosen with the range and the spacing reported in Table 1 . The set of parameters characterised by h MS min ≥ 1 were excluded since this choice does not produce an action potential. The statistics of the pacemaker behaviour was thus evaluated on a total of 52598 samples. For both models, an external stimulus was applied at ) and ( ) represent the points 1,2 and 3 for the mMS and MS models, respectively. The zeroth point ( •) represents the rest (initial) state and coincides for both models. The values of the parameters were taken from [16] . t = 0 ms and the numerical solution was evaluated over 1200 ms; parameter sets presenting more than one depolarisation were considered pacemaker. For 3042 combinations ( 6% of the cases) MS displayed pacemaker behaviour, while the mMS was robust in all of the 52,598 tests.

Asymptotic derivation of the restitution curve
Following the same procedure adopted in [16] , it is possible to derive an explicit leading order asymptotic approximation for the APD restitution curve of mMS, based on the assumption τ in τ out τ open , τ close (9) From assumption (9) it follows that the time constants of Eq. (3) are infinitesimal if compared to the time constants of Eq. (4) and thus the duration of the depolarisation and repolarisation phases can be neglected when compared to the APD and recovery durations.

Single stimulus: maximum APD (APD max )
In this paper we adopt the same definition introduced in [16] for APD as the elapsed time during which v m ≥ v gate .
Due to the separation of the time scale introduced by assumption (9) and considering a cell membrane initially in the rest state

Multiple stimuli: restitution curve
In this section we derive a leading order approximation of the APD when two (or more) stimuli are applied. The first stimulus is applied at t = 0 ms to a fully recovered tissue; the second stimulus is applied at time t = S ms , to a system not fully recovered.
Denoting by APD 1 the APD following the first stimulus and by DI = S-APD 1 the diastolic interval (DI), that is the time elapsed between the end of the APD and the second pacing stimulus, from Eq. (4) there follows Due to assumption (9) , the depolarisation duration is negligible if compared to APD and recovery duration; thus, from (4) and the initial condition (11) , the following relation between APD and DI, called restitution, holds: lines depict the APD restitution computed by numerically solving MS and mMS. In Fig. 3 B, the same restitutions are plotted for v gate = v * gate = 1 − 1 − h MS min : this value is evaluated by imposing h mMS min = h MS min and produces identical analytical restitution curves for both models. This choice of the activation threshold also produced similar computed restitution curves for the two models (dashed lines, Fig. 3 B).
The restitution curves were computed by applying an s 1 _ s 2 protocol [17] . Briefly, the model is repetitively stimulated with an inter-pacing interval s 1 until a steady state is achieved, (in the present paper, s 1 = 10 0 0 ms , applied 10 0 times to achieve a limit cycle), and then a single stimulus is applied at an interval s 2 < s 1 and the system is left to evolve. This procedure is applied for different values of s 2 and the APD is measured after the premature stimulus is applied. In the present paper, the values of the premature stimulus s 2 < s 1 were chosen with the following sequence: (10) and (12) are formally equivalent to those obtained in [16] for MS and differ only in the expression of Point 2 defined in Section 3.1 . While in the MS model this point is defined by (5) and does not depend on v gate , in the mMS model this point is defined by expression (7) which introduces a dependence on v gate in expressions (10) and (12) .
Since the APD is defined as the elapsed time during which v m ≥ v gate , expressions (10) and (12) are able to account for this dependency. As a demonstration, if v gate = 1 . 0 is chosen, one expects no APD will be generated, independently of the choice of the other 4 ionic parameters. The analytic restitution expression of the MS model still will furnish a non-zero APD in this case. In contrast, from expressions (10) and (11) the restitution equation for the mMS model predicts an APD of 0 ms.

Numerical examples
To compare the functional characteristics of the MS and mMS models we compare model simulations for a stable action potential, an action potential affected by pacemaker behaviour, a dynamic restitution protocol and pacemaker behaviour in a tissue simulation.
Time discretisation was performed by a backward Euler method and non-linearities were treated by Newton iterations.
A time step dt = 0 . 005 ms was chosen; to test the accuracy of the chosen time step we evaluated the numerical solution for the set of parameters reported in [16] with a time step dt = 0 . 005 ms and a time step dt = 0 . 0 0 05 ms . The maximum and the L 2 differences between the solutions were equal to 0.004 and 0.005 respectively, (i.e., 0.5% of the maximum value of v m ). The externally applied stimulus is characterised by an intensity of 1 . 0 ms −1 and a duration of 0.4 ms.
In Section 5.1 we deal with the ionic parameters taken from [16] that do not provide any pacemaker activity for MS.
In Section 5.2 we adopt the parameters reported in Table 2 that yield pacemaker behaviour for the MS model. We will show that while the MS model will furnish a pacemaker solution, the mMS model will be robust to this behaviour. In Section 5.3 we evaluate the dynamic restitution curves for both models with the parameter set defined in [16] and we show that the differences between the two models in the dynamic APD restitutions are small.
In Section 5.4 we consider a homogeneous tissue slab characterised by the parameters reported in Table 3 , stimulated by a cross-field protocol. This set of parameters did not produce any pacemaker behaviour for the MS model in a single cell, nor in a 1D tissue string for the same parameter set and model conditions, but did show pacemaker behaviour in two dimensions. We show that the mMS model is robust to pacemaker behaviour under these conditions.

Example 1: pacemaker free MS
In the first example the MS and the mMS models were solved with the parameter set taken from [16] . This choice yields a pacemaker free MS model. Three stimuli were applied at an inter-pacing interval T = 700 ms . For both models the transmembrane potential ( Fig. 4 A) and the gate variable ( Fig. 4 B) functions of time are depicted.
The dynamics of the gate variable h do not significantly differ between the two models, while the trans-membrane potential v m obtained by solving the mMS model only differs in the maximum value reached at the end of the depolarisation.
In Fig. 5 the trans-membrane potential v m and the gate variable h are plotted when a value of v gate = v * gate is chosen. This choice of v gate yields the same analytical restitution curves for the mMS and the MS models.

Example 2: MS with pacemaker activity
In the second example we characterise both models with the parameter set summarised in Table 2 . This choice leads to a pacemaker solution in the MS model. The system is paced once at t = 0 and the solution is evaluated until t = 1200 ms . The results are plotted in Fig. 6 ; while the mMS model presents only one     Table 2 .  Table 2 were taken.
action potential, the MS cyclically depolarises (with 4 activations during the simulation period). In Fig. 7 the nullclines and the phase portrait are depicted for both models. The mMS model ( Fig. 7 A) recovers without crossing the left nullcline, since the recovery of h is far from the nullcline. The MS model ( Fig. 7 B) crosses the left nullcline during recovery and thus re-depolarises entering into a periodic limit cycle.

Example 3: dynamic restitution curves
Dynamic restitution curves are used to characterise the value of the APD when the tissue is periodically paced with a pacing period S . As discussed in [16] , at larger values of S , each stimulus produces an action potential, yielding a 1:1 correspondence between S and APD, called 1:1 behaviour. As S decreases, this 1:1 behaviour eventually becomes unstable; as a result, two possible behaviours can be presented: • A 2:1 behaviour: only every other stimulus produces an action potential; • A 2:2 behaviour (alternans): each stimulus will produce an action potential, but the APD periodically alternates between short and long.
For both models, dynamic restitution curves were numerically evaluated. For each pacing period S , the system was periodically paced by applying 104 external stimuli and the APD was computed for each of the last 4 applied stimuli. When the APD values of the last 4 stimuli all coincide, the behaviour was considered of type 1:1. If the APD presented two different alternating values and all of the 4 stimuli produced an action potential the behaviour was considered of type 2:2 (alternans). Lastly, if the APD presented only one value and the 4 stimuli produced only two activations, the behaviour was considered to be of type 2:1. The examples tested did not present any type 2:1 behaviour.
Dynamic restitutions were evaluated for the parameters of Table 2 and for v gate = v * gate as defined in Section 4.2 . The pacing interval S was decremented from 700 ms to 300 ms with a step of 100 ms and then from 280 ms with a step of 2 ms to the first value that did not produce APD. For v gate = 0 . values differing by 240 ms in the MS model. The pacing periods where either APD bifurcates or no APD was produced differ by 2 ms between both models, a difference comparable with the step used to decrement S .
In Fig. 8 A the dynamic restitution of the MS model is depicted, while in Fig. 8 B the same curve is depicted for the mMS model. In Fig. 8 C and D the dynamic restitutions are depicted for the mMS and the MS models respectively. The two different APD values are depicted with a blue circle and a red triangle (second stimulus). When the circles and triangles are no longer superimposed, bifurcation occurs.

Example 4: 2D model with a cross-field stimulus
In this example we consider an homogeneous isotropic tissue slab measuring 5 cm × 5 cm; the tissue electrophysiology is described by the mono-domain model [13] and solved with the CARP [18,26] finite element code. The equations were discretised in space with a triangular mesh with a characteristic size h = 189 μm and in time with a time step of 0.1 ms. The model was simulated for t = 3500 ms .
The tissue electrophysiology was characterised by the parameters reported in Table 3 . These parameters did not generate any pacemaker activity in simulations of an isolated single cell or in a 1D domain.
The tissue was stimulated by the cross field stimulation protocol, [19] . Briefly, an external stimulus was applied to two orthogonal regions at two different stimulation times to trigger a spiral activation pattern. In this work the first stimulus s 1 was applied at t = 0 ms on the left edge of the slab between x = 0 and x = 0.2 cm, and s 2 was applied at t = 290 ms on the bottom edge of the slab between y = 0 and y = 0.25 cm. For both stimuli, a current intensity of J stim = 2 . 0 ms −1 with a duration of T stim = 0.6 ms were employed.
The transmembrane potentials of the mMS (top row) and the MS (bottom row) models are depicted in Fig. 9 for t = 290 ms ( s 2 is applied), t = 390 ms (a spiral wave is generated, no pacemaker activity is present), t = 420 ms (pacemaker behaviour appears in the MS model) and t = 620 ms (presence of pacemaker behaviour in the MS model). A video of the whole simulation can be found in the online supplement.
The mMS model generated a spiral re-entry wave that broke up at t = 2080 ms and terminated at t = 2470 ms , whereas the MS model showed a pacemaker behaviour, initiated by a pacemaker beat at t = 420 ms ( Fig. 9 , third column) and which persisted over the entire simulation period.  . 9. Trans-membrane potential followed by a cross-field stimulus for the mMS (top) and the MS (bottom) models at times t = 290, 390, 420 and 620 ms.

Discussion and conclusions
In this paper we introduced a novel two state-variable ionic model, obtained by modifying the standard Mitchell-Schaeffer (MS) ionic model. The new model has the same benefits as the standard MS model and is proven to be robust to pacemaker behaviour. Previously the MS model had been adapted to introduce pacemaker behaviour [6] . The nonlinear term characterising the inward ionic current was modified by introducing two new parameters. In this paper, we presented a model that is robust to pacemaker behaviour by modifying the inward current in a similar manner. Unlike [6] , in this paper we also gated the outward current. This yields three advantages: first, the threshold value of the transmembrane voltage above which an action potential is triggered corresponds to v gate ; second, at the end of the depolarisation, v m = 1 ; third, the analytical solution of the mMS model coincides with the analytical solution of the MS model for the same set of ionic parameters and v gate = v * gate . We then introduced an asymptotic derivation of the restitution curves, obtaining a relation formally equivalent to that obtained in [16] . We also showed, through numerical examples, that there are only small differences in the action potential introduced by the modified MS (mMS) equations compared to the original MS model when the same set of parameters are adopted. We also compared the APD dynamic restitution curves obtained by numerically solving both models; we showed that the mMS is able to reproduce the APD alternans and that these occur under a similar pacing regime in both models. Last, we demonstrated the robustness of the new model even when incorporated into a mono-domain tissue simulation, confirming its applicability in tissue scale patient-specific modelling.
Research Centre award to Guy's and St Thomas' NHS Foundation Trust in partnership with King's College London and King's College Hospital NHS Foundation Trust.

Supplementary material
Supplementary material associated with this article can be found, in the online version, at 10.1016/j.mbs.2016.08.010