Royal Society Publishing


We analyse small parameters in selected models of biological excitability, including Hodgkin–Huxley (Hodgkin & Huxley 1952 J. Physiol. 117, 500–544) model of nerve axon, Noble (Noble 1962 J. Physiol. 160, 317–352) model of heart Purkinje fibres and Courtemanche et al. (Courtemanche et al. 1998 Am. J. Physiol. 275, H301–H321) model of human atrial cells. Some of the small parameters are responsible for differences in the characteristic time-scales of dynamic variables, as in the traditional singular perturbation approaches. Others appear in a way which makes the standard approaches inapplicable. We apply this analysis to study the behaviour of fronts of excitation waves in spatially extended cardiac models. Suppressing the excitability of the tissue leads to a decrease in the propagation speed, but only to a certain limit; further suppression blocks active propagation and leads to a passive diffusive spread of voltage. Such a dissipation may happen if a front propagates into a tissue recovering after a previous wave, e.g. re-entry. A dissipated front does not recover even when the excitability restores. This has no analogy in FitzHugh–Nagumo model and its variants, where fronts can stop and then start again. In two spatial dimensions, dissipation accounts for breakups and self-termination of re-entrant waves in excitable media with Courtemanche et al. kinetics.


1. Introduction

The motivation of this study comes from a series of numerical simulations of spiral waves (Biktasheva et al. 2003) in a model of human atrial tissue based on the excitation kinetics of Courtemanche, Ramirez & Nattel 1998 (CRN). The spiral waves in this model tend to break up into pieces and even spontaneously self-terminate (see figure 1).

Figure 1

Self-termination of a spiral wave in the CRN model. Dark-shaded component: higher values of the transmembrane voltage E, i.e. excitation; light-shaded component: gating variable Embedded Image. Diffusion coefficient Embedded Image, preparation size 75×75 mm. See also Biktasheva et al. (2003).

No mathematical model of cardiac tissue is now considered ultimate or can claim absolute predictive power. The spontaneous self-termination may be relevant to human atrial tissue or may be an artefact of modelling. Understanding the mechanism of this behaviour in some simple terms would allow a more direct and certain verification. This is difficult as the models are very complex and the events depicted in figure 1 have many different aspects. Traditionally, such understanding has been achieved in terms of simplified models, starting from axiomatic cellular automata description (Wiener & Rosenblueth 1946) to simplified partial differential equation (PDE) models (FitzHugh 1961; Nagumo et al. 1962), which allow asymptotic study by means of singular perturbation techniques (Tyson & Keener 1988). This approach can describe some of the features observed in figure 1, e.g. the ‘action potential duration restitution slope 1’ theory predicts when the stationary rotation of a spiral wave is unstable against alternans (Nolasco & Dahlen 1968; Karma et al. 1994; Karma 1994). The relevance of the ‘slope-1’ theory to particular models is debated (Cherry & Fenton 2004), but in any case it only predicts the instability of a spiral wave, not whether it will lead to complete self-termination of the spiral wave, its breakup, or just meandering of its tip. We need to understand how the propagation of a wave is blocked. This has unexpectedly turned out to be rather interesting. Some features of the propagation block in figure 1 can never be explained within the standard FitzHugh–Nagumo (FHN) approach. As this was the only well-developed asymptotic approach to excitable systems around, we had either to accept that this problem is too complicated to be understood in simplified terms, or to develop an alternative type of simplified model and corresponding asymptotics.

We chose the latter. This paper summarizes our progress in this direction in the last few years (Biktashev 2002, 2003; Biktasheva et al. 2003; Suckley & Biktashev 2003; Biktashev & Suckley 2004; Suckley 2004; Biktashev & Biktasheva 2005). The results on the asymptotic structure of the CRN model are published for the first time.

2. Tikhonov asymptotics

The standard Tikhonov–Pontryagin singular perturbation theory (Tikhonov 1952; Pontryagin 1957) summarized in Arnold et al. (1994) is usually formulated in terms of ‘fast–slow’ systems in one of the two equivalent forms,Embedded Image(2.1)where ϵ> 0 is small, t is the ‘slow time’, T is the ‘fast time’, Embedded Image, X is the vector of ‘slow variables’ and Y is the vector of ‘fast variables’. The theory is applicable if all relevant attractors of the fast subsystem are asymptotically stable fixed points. So the variables have to be explicitly classified as fast and slow, and the system should contain a small parameter which formally tends to zero although the original system is formulated for a particular value of that parameter, say 1.

We consider Hodgkin & Huxley (1952) (henceforth referred to as HH) and Noble (1962) (henceforth N62) models. Both can be written in the same form asEmbedded Image(2.2)where Embedded Image is the total transmembrane current (μA cm−2), t is time (ms), E is the transmembrane voltage (mV), Embedded Image, Embedded Image are the reversal potentials of sodium, potassium and leakage currents, respectively (mV), Embedded Image are the corresponding maximal specific conductances (mS cm−2), n, m, h are dimensionless ‘gating’ variables, Embedded Image is the specific membrane capacitance (μF cm−2), Embedded Image are the gates' instantaneous equilibrium ‘quasi-stationary’ values and Embedded Image are the characteristic time-scale coefficients of the gates dynamics (ms). Definition and comparison of parameters and functions used in (2.2) for HH and N62 models can be found in Suckley & Biktashev (2003).

To classify the dynamic variables by their speeds, we define empiric characteristic time-scale coefficients, Embedded Image. For a system of differential equations Embedded Image, Embedded Image, we define Embedded Image. The τ's obtained for m, h and n in this way coincide with Embedded Image in (2.2), and this definition can be extended to other variables, e.g. E in the case of (2.2).

(a) Hodgkin–Huxley

Figure 2a shows how τ's change during a typical action potential (AP) in the HH model. The speeds of E and m exchange places during the AP, as do the speeds of h and n, but at all times the pair Embedded Image remains faster than the pair Embedded Image. This suggests the introduction into system (2.2) of a parameter ϵ, which in the limit Embedded Image makes variables E and m much faster than n and h:Embedded Image(2.3)

Figure 2

Tikhonov asymptotics of HH. (a) Characteristic times Embedded Image during an AP. Thin solid magenta line: the AP for a reference. (b) Phase portrait of the fast subsystem (2.4) at Embedded Image, Embedded Image. Solid red line: horizontal isocline; dashed blue line: vertical isocline; filled black circles: stable equilibria; dash-dotted green line: stable separatrix of the saddle, the boundary of attraction basins; black dotted lines: selected trajectories. (c) A three-dimensional view of the slow manifold (the surface). Solid line: the fold line; dotted line: the selected trajectory and its projections on coordinate walls. (d) AP in the original model, ϵ=1, solid red line, and when the fast variables are made faster, Embedded Image, dashed blue line. See also Suckley & Biktashev (2003).

The properties of this system in the limit Embedded Image are shown in figure 2b,c. The fast transient, corresponding to the AP upstroke, can be studied by changing the independent time variable in (2.3) to Embedded Image and then considering the limit Embedded Image, which gives the fast subsystem of two equations for m and E,Embedded Image(2.4)in which the slow variables h and n are parameters as their variations during the onset are negligible. An example of a phase portrait of system (2.4) at selected values of h and n is shown in figure 2b. It is bistable, i.e it has two asymptotically stable equilibria, and a particular trajectory approaches one or the other depending on the initial conditions. The basins of attraction of the two equilibria are separated by the stable separatrices of a saddle point, which is the threshold between ‘all’ and ‘none’ responses. A fine adjustment of initial conditions at the threshold will cause the system to come to the saddle point. This is a mathematical representation of the excitation threshold in Tikhonov asymptotics.

For different values of n and h, the location of the equilibria in the fast subsystem vary. All equilibria Embedded Image at all values of n and h form a two-dimensional slow manifold in the four-dimensional phase space of (2.3) with coordinates Embedded Image. Projection of this two-dimensional manifold into the three-dimensional space with coordinates Embedded Image is depicted in figure 2c. It has a characteristic cubic folded shape, with two fragments of a positive slope (as it appears on the figure), separated by an ‘overhanging’ fragment of a negative slope. The borders between the fragments are the fold lines, seen as nearly horizontal solid curves on the picture. The positive slope fragments consist of stable equilibria, and the negative slope fragment consists of unstable equilibria (saddle points) of the fast subsystem.

The points of this manifold are steady-states if considered on the time-scale Embedded Image or equivalently Embedded Image. On the time-scale Embedded Image, these points are no longer steady states, but we observe a slow (compared to the initial transient) movement along this manifold, which explains its name. Asymptotically, the evolution on the scale Embedded Image can be described by the limit Embedded Image in (2.3), which gives a system of two finite equations and two differential equations,Embedded Image(2.5)The finite equations define the slow manifold and the differential equations define the movement along it.

Figure 2c shows a selected trajectory of system (2.3) corresponding to a typical AP solution. The only equilibrium of the full system (2.3), corresponding to the resting state, is at the lower, ‘diastolic’ branch of the slow manifold. If the initial condition is in the basin of the upper branch, the trajectory starts with a fast initial transient, corresponding to the upstroke of the AP, then proceeds along the upper, ‘systolic’ branch of the slow manifold, which corresponds to the plateau of the AP. When the trajectory reaches the fold line, a boundary of the systolic branch, the plateau stage is over. The fast subsystem is no longer bistable. The only stable equilibrium is now at the diastolic branch. So the trajectory jumps down. This is repolarization from the AP and it happens at the time-scale Embedded Image. The trajectory then slowly proceeds along the diastolic branch towards the resting state.

So, an inevitable feature of the asymptotics (2.3) is that the solution at smaller ϵ has not only a faster upstroke, but also a faster repolarization, and the asymptotic limit of the AP shape is rectangular as opposed to the triangular shape in the exact model, see figure 2d. This is undesirable as it means that asymptotic formulae obtained in this way produce qualitatively inappropriate results.

The practical importance of this exercise is limited, as in HH the AP are not much longer than upstrokes.

(b) Noble (1962)

This model is more relevant to cardiac AP. The speed analysis of N62 model, similar to the one we have done for HH model, reveals a different asymptotic structure but ultimately similar results. Figure 3a demonstrates three rather than two different time-scales. Variable m is the fastest of all, we call it ‘superfast’. Of the remaining three, variables E and h are fast and variable n is slow. So we need two small parameters, Embedded Image to describe the difference between the fast and superfast time-scales, and Embedded Image to describe the difference between the slow and fast time-scales. System (2.2) then takes the formEmbedded Image(2.6)Consider first the limit Embedded Image. The superfast subsystem consists of one differential equation for m. It always has exactly one equilibrium which is always stable. So after a supershort transient, Embedded Image is always close to its quasi-stationary value Embedded Image. Thus, with an error of the order of Embedded Image, we may approximate m by Embedded Image and discard the first equation, i.e. adiabatically eliminate superfast variable m.

Figure 3

Tikhonov asymptotics of N62. (a, b, d) Notations similar to figure 2. (c) A two-dimensional view of the stable (black solid) and unstable (green dashed) branches of the slow manifold, and typical pacemaker potential trajectories (the limit cycles): solid red (Embedded Image) and dashed blue (Embedded Image), vertical dash-dotted line shows value Embedded Image for the fast phase portrait on (b). See also Suckley & Biktashev (2003).

With the remaining system of three differential equations for E, h and n, we consider the change of the time variable as before, Embedded Image, and proceed to the limit Embedded Image. This produces the fast subsystem in the formEmbedded Image(2.7)Figure 3b shows a phase portrait of this system for a selected value of n when (2.7) is bistable. As there is one slow variable, all the equilibria of the fast subsystem form a one-dimensional manifold, i.e. a curve, in the three-dimensional phase space with coordinates Embedded Image. Its projection on the plane Embedded Image is shown in figure 3c. Again the stable equilibria correspond to one slope (negative for the given choice of coordinates) and the opposite slope corresponds to the unstable equilibria. The branches are separated from each other by fold points. Again we have an upper, systolic branch, separated from the lower, diastolic branch, and the system has no alternative but to jump from one branch to the other in the time-scale Embedded Image. As it happens, there are no true equilibria in the N62 model at standard parameter values, so these jumps happen periodically, producing pacemaker potential. This corresponds to the automaticity of cardiac Purkinje cells. Certain physiologically feasible changes of parameters may produce an asymptotically stable equilibrium at the diastolic branch, i.e. turn an automatic Purkinje cell into an excitable cell.

The systolic branch is separated from the diastolic branch, so in the asymptotic limit Embedded Image, if the upstrokes are fast, the repolarizations are similarly fast. This is in contradiction with the behaviour of the full model, which makes such an asymptotic analysis unsuitable.

3. Non-Tikhonov asymptotics

(a) Noble (1962)

To overcome the difficulties of the Tikhonov approach, we have developed an alternative based on actual biophysics behind N62 and other cardiac excitation models, see figure 4a. The upstroke of an AP is much faster than the repolarization, as the two processes are caused by different ionic currents. The upstroke is very fast as the fast Na current causing it is very large. However, all other currents, including the outward K currents that bring about repolarization, are not as large, and there is no reason to tie the two classes of currents together in an asymptotic description. Mathematically, this means that in the right-hand side of the equation for E, only the term corresponding to the fast Na current is large and should have the coefficient Embedded Image in front of it.

Figure 4

Non-Tikhonov asymptotics of N62, excitable variant. (a) Functions of E illustrating the small quantities taken into account. (b) Phase portrait of the fast subsystem. (c) Phase portrait of the slow subsystem. (d) Action potential solution of (3.2) for Embedded Image, solid red, and for Embedded Image, dashed blue. See also Biktashev & Suckley (2004).

Next, the fast Na current is large only during AP upstroke, and remains small during other stages. The current is regulated by two gating variables, m and h, and the quasi-stationary value of the specific conductivity of the current, defined as Embedded Image, is always much smaller than 1. This happens because Embedded Image and Embedded Image is very small for E below a certain threshold voltage Embedded Image, and Embedded Image is very small for E above another threshold voltage Embedded Image and Embedded Image. So the ranges of almost complete closure of Embedded Image and Embedded Image overlap. So whenever E changes so slowly that m and h have enough time to approach their quasi-stationary values, the fast Na channels are mostly closed. The possibility for the opening of a large fraction of Na channels only exists during the fast upstroke, as m gates are much faster than h gates and have time to open before h closes.

Thus, the facts that m gate is much faster than h gate, Embedded Image for Embedded Image, Embedded Image for Embedded Image and Embedded Image, are all related and reveal why the upstroke of the AP is much faster than all other stages.

We adopt the hierarchy of times suggested by the formal speed analysis in §2, i.e. m is a superfast variable, E and h are equally fast variables and n is a slow variable. We keep the same notation Embedded Image and Embedded Image for the corresponding small parameters. The small parameters should also ensure that Embedded Image and Embedded Image are small in some ranges of E but not in others. We identify this smallness with Embedded Image rather than Embedded Image, as it is to compensate the large value of Embedded Image outside the upstroke, and the large value of Embedded Image is described by Embedded Image. We denote the Embedded Image-dependent versions of Embedded Image and Embedded Image as Embedded Image and Embedded Image. In this way, we obtainEmbedded Image(3.1)where Embedded Image is the Heaviside function.

The limit Embedded Image gives adiabatic elimination of the m gate, Embedded Image.

The analysis of the limit Embedded Image is complicated by a feature incidental to N62 and not found in other cardiac excitation models. The small conductivity of the window current, Embedded Image, is multiplied by a large factor Embedded Image. In N62, the resulting window component of Embedded Image is comparable to other small currents and cannot be neglected outside the upstroke. The implications of this complication are analysed in Biktashev & Suckley (2004). The result, in brief, is that the Embedded Image-dependent part of (3.1) can, in the limit Embedded Image, be replaced with a ‘modified N62’ model:Embedded Image(3.2)where Embedded Image, and as before Embedded Image for Embedded Image, Embedded Image for Embedded Image and Embedded Image.

The limit Embedded Image of (3.2) in the fast time Embedded Image gives the fast subsystemEmbedded Image(3.3)As intended, (3.3) takes into account only the fast sodium current and the gates controlling it, and everything else is a small perturbation on this time-scale. The phase portrait of (3.3) is unusual, see figure 4b. The horizontal isocline (the red set) is not just a curve but contains a whole domain Embedded Image. The vertical isocline (the blue set) lies entirely within the red set, so the whole line Embedded Image consists of equilibria. An upstroke trajectory may end up in any of the equilibria above Embedded Image, so the height of an upstroke depends on initial conditions. For subthreshold initial condition, voltage remains unchanged in the fast time-scale. Exactly what happens at the threshold Embedded Image depends on details of approximating function Embedded Image, but in any case it does not involve any unstable equilibria. This is all different from Tikhonov systems (see the paragraph after equation (2.4)) where the height of the upstroke is fixed, subthreshold potential decays in the fast time-scale and the threshold consists of unstable equilibria and, if appropriate, their stable manifolds. In this sense, asymptotics of (3.3) give a new meaning to the notion of excitability, completely different from that in the Tikhonov systems.

Let us consider the slow subsystem of (3.2). For any value of n, we have a whole line of equilibria in the fast system Embedded Image. The collection of such lines makes a two-dimensional manifold in the three-dimensional space with coordinates Embedded Image. So the fast variable h can be adiabatically eliminated on the time-scale Embedded Image. Thus, the slow subsystem, i.e. the limit Embedded Image in (3.2), isEmbedded Image(3.4)The phase portrait of this system is shown in figure 4c. Further discussion of its properties can be found in Biktashev & Suckley (2004). Note that voltage E features in both the fast and slow subsystems, i.e. it is a fast or a slow variable depending on circumstances. This kind of behaviour is not allowed in Tikhonov asymptotic theory, so it is ‘a non-Tikhonov’ variable.

(b) Courtemanche et al. (1998)

CRN is a system of 21 ordinary differential equations modelling electric excitation of human atrial cells (see Courtemanche et al. (1998) for a description).

Formal analysis of the time-scales Embedded Image of dynamic variables by the same method as we used for HH and N62, reveals a complicated hierarchy of speeds, which changes during the course of the AP (see figure 5a). From variables with smaller τ's, we select those that remain close to their quasi-stationary values during an AP, and which can be replaced by those quasi-stationary values without significantly affecting the AP solution. We call them superfast variables. These include m, Embedded Image and w. As before, we denote the associated small parameter Embedded Image.

Figure 5

Non-Tikhonov asymptotics of the Courtemanche et al. (1998) model. (a) Characteristic times during a typical AP potential solution. Red solid line: Embedded Image; blue dashed line: τ's of superfast (m, Embedded Image, w) and fast (h, Embedded Image, d); green dash-dotted lines: τ's of other, slow variables. (b) Phase portrait of the fast subsystem (3.6). (c) Action potential upstroke, Embedded Image (solid red), Embedded Image, Embedded Image (dotted blue), Embedded Image (dashed green), in the fast time Embedded Image. (d) Same, for the whole AP in the slow time t. See also Suckley (2004).

Next, we identify the fast variables with speeds comparable to the AP upstroke. This is also done by comparing the instantaneous values of the variables with the corresponding quasi-stationary values, and checking how their adiabatic elimination affects the AP, for the AP solution after the initial upstroke. In this way, we identify variables h, Embedded Image and d as fast, with associated small parameter Embedded Image.

Similar to N62, the transmembrane voltage is Embedded Image-fast only during the AP upstroke due to Embedded Image-large values of Embedded Image during that period, and is slow otherwise. This is due to nearly perfect switch behaviour of the gates m and h. The definition of Embedded Image in this model is more complicated as there is also the j gate; however, j is slow and does not change noticeably during the upstroke.

These considerations lead to the following parameterization of the model:Embedded Image(3.5)where Embedded Image is the fast Na current and Embedded Image represents all other currents. Here we have shown only equations that contain Embedded Image or Embedded Image. All other equations are the same as in the original model.

As before, we adiabatically eliminate the superfast variables in the limit Embedded Image, and turn Embedded Image in the fast time-scale Embedded Image. This gives the fast subsystemEmbedded Image(3.6)Note that in (3.6) the equations for Embedded Image and d depend on E, but equations for E and h do not depend either on Embedded Image or on d. Thus, the evolution equations for E and h form a closed subsystem. This system is identical to (3.3), up to the values of parameters, definitions of the functions of E and the presence of the slow variable j as the factor at the maximal conductance of the Na current, Embedded Image. The phase portrait of (3.6), figure 5b, is similar to that of N62, figure 4b. So the peculiar features of the fast subsystem of N62 are not unique and are found in many cardiac models, including CRN.

With a view of a practical application of approximation (3.6), it is interesting to test its quantitative accuracy. This is illustrated in figure 5c for the shape of the upstroke in the fast time Embedded Image and figure 5d for the shape of the AP in the slow time t. We see that the approximation of the AP is very good in both limits Embedded Image and Embedded Image, except for the upstroke: e.g. the peak voltage is overestimated by about 13 mV. This is mostly due to the limit Embedded Image, i.e. replacement of m with Embedded Image. So the accuracy could be significantly and easily improved by retracting the limit Embedded Image, which amounts to inclusion of the evolution equation for m instead of the finite equation Embedded Image. Most of the qualitative analysis remains valid. However, here for simplicity, we stick to the less accurate but simpler case Embedded Image. Note that in figure 5b–d, we took Embedded Image, Embedded Image; the error introduced by that was small compared to other errors, particularly the error introduced by Embedded Image.

4. Application to spatially distributed systems

(a) Fronts

We now use the fast Na subsystem of the cardiac excitation (3.6) to consider a propagation of an excitation front through a cardiac fibre. In one spatial dimension, this requires replacement of ordinary time derivatives with partial derivatives and adding a diffusion term into the equation for E:Embedded Image(4.1)where Embedded Image and we have put Embedded Image. We consider solutions in the form of propagating fronts. For definiteness, let us assume the fronts propagating leftwards, so Embedded Image, Embedded Image, Embedded Image, Embedded Image, Embedded Image, Embedded Image and Embedded Image. In this formulation, we have three constants characterizing the solutions, the pre-front voltage Embedded Image, the post-front voltage Embedded Image and the speed c. It is not obvious which combinations of the three parameters admit how many front solutions. So we have considered a ‘caricature’ of (4.1) by replacing functions Embedded Image and Embedded Image in it with constantsEmbedded Image(4.2)

This system is piecewise linear and admits complete analytical investigation. Details can be found in Biktashev (2002, 2003); here we only briefly outline the results. Figure 6a illustrates a typical front solution. It exists if speed c and pre-front voltage Embedded Image satisfy a finite equation involving also the constants J and Embedded Image. The resulting dependence of the conduction velocity c on excitability J for a few selected values of Embedded Image is shown in figure 6b. These front solutions exist only for J at or above a certain minimum Embedded Image which depends on Embedded Image. For Embedded Image, there are two solutions with different speeds. Numerical simulations of PDE system (4.2) suggest that solutions with higher speeds are stable and solutions with lower speeds are unstable; this has been confirmed analytically by Hinch (2004).

Figure 6

Fronts in the spatially distributed Na subsystem. (a) The structure of the front solution in the caricature model. (b) Speed of the fronts as a function of excitability, at selected values of the pre-front voltage, in the caricature model (4.2). Embedded Image is minimal excitability at which propagation is possible at any Embedded Image, and Embedded Image is the corresponding propagation speed. (c) Speed of the fronts as a function of excitability, at selected values of pre-front voltage, in the Na subsystem of the CRN model (4.1). On (b) and (c), solid red lines, above and rising, are the stable branches and dashed green lines, below and decreasing, are the unstable branches. (d) For comparison: speed of the fronts in a typical Tikhonov front (fast susbsystem of the FHN model). See also Biktashev (2002) and Biktashev & Biktasheva (2005).

The replacement of functions Embedded Image and Embedded Image with constants is a rather crude step. The purpose of the caricature is not to provide a good approximation, but to investigate qualitatively the structure of the solution set. To see if this structure is the same for the more realistic models, we have solved numerically the boundary-value problems for the front solutions in (4.1). There the role of the excitability parameter is played by the variable j. The results of the calculations are shown in figure 6c. Not only is the topology of the solution set the same, but the overall behaviour of Embedded Image in (4.1) is quite similar to that of Embedded Image in (4.2), despite the crudeness of the caricature.

PDE simulations show that approximation (4.1) overestimates the conduction velocity by almost 50% compared to the full model, and the error is again mainly due to the adiabatic elimination of the m gate.

After eliminating the superfast variables m, Embedded Image and w and the fast variables h, Embedded Image and d, and retaining the non-Tikhonov variable E, the slow subsystem of (3.5) has 16 equations. It describes the AP behind the front.

The most important conclusion is that for any particular value of the pre-front voltage Embedded Image, there is a certain minimum excitability Embedded Image and corresponding minimum propagation speed Embedded Image, and for Embedded Image, no steady front solutions are possible. This is completely different from the behaviour in FHN type systems, where local kinetics are Tikhonov and a front is a trigger wave in a bistable reaction–diffusion system. A typical dependence of the speed of such a trigger wave on a slow variable is shown in figure 6d: it can be slowed down to a halt or even reversed. The reversed trigger waves describe backs of propagating pulses in FHN systems. Thus, questions about the shape of the backs of APs and propagating pulses and the spectrum of propagation speeds of a propagating excitation wave in a tissue come to be closely related. In both questions, our new non-Tikhonov approach provides different answers from the traditional Tikhonov/FHN approach. We have already seen that the new description is more in line with the detailed ionic models regarding the back of an AP. In §4b, we demonstrate the advantage of the new approach regarding the fronts.

(b) Dissipation of fronts

The fast subsystem of a typical spatially dependent cardiac excitation model, discussed in §4a, only provides part of the answer. This description should be completed with the description of the slow movement. The fronts are passing so quickly through every given point that the values of the slow variables at that point change little while it happens. Away from the fronts, the fast variables keep close to their quasi-stationary values. In our asymptotics, this means, in particular, that the fast Na channels are closed, and E is not a fast but a slow variable. Assuming absence of spatially sharp inhomogeneities of tissue properties, simple estimates show that outside fronts, the diffusive current is much smaller than ionic currents, so dynamics of cells outside fronts are essentially the same as dynamics of isolated cells outside AP upstrokes.

Propagation of the next front depends on the transmembrane voltage E, which serves as parameter Embedded Image, and the slow inactivation gate j of the fast Na current. This dependence gives an equation of motion for the front coordinate Embedded Image,Embedded Image(4.3)where the instantaneous speed of the front c is determined by the values of E and j at the sites through which the front traverses (in the case of E this is the value which would be there if not the front). This can only continue as long as the function Embedded Image remains defined, i.e. while Embedded Image. If the front runs into a region where this is not satisfied, its propagation becomes unsustainable.

What will happen then is illustrated in figure 7a, where the parameter Embedded Image was varied in space and time. To make the effect more prominent, we did not use smooth variation, but put Embedded Image in the left half of the interval for some time and then restored it to its normal value. The propagating front reached this region while it was in the inexcitable state. The result was that the sharp front ceased to exist, it ‘dissipated’, and instead of an active front we observed a purely diffusive spread of the voltage. The excitability was restored a few milliseconds later, but the sharp front did not recover and diffusive spread of voltage continued, leading eventually to a complete decay of the wave. Note that the back of the propagating pulse was still very far when the impact that caused the front dissipation happened.

Figure 7

Effects of temporary propagation blocks. (a) In the CRN model. (b) In the FHN model. (c) In the caricature Na front. (d) In the Na subsystem of the CRN model. See also Biktashev & Biktasheva (2005).

This is completely different from the behaviour of a FHN system in similar circumstances, shown in figure 7b. There, propagation was blocked for almost the whole duration of the AP. And yet when the block was removed, the propagation of the excitation wave resumed. Only if the block stays so long that the waveback reaches the block site and the ‘wavelength’ reduces to zero, the wave would not resume. Such considerations have led to a widespread, would-be obvious assumption that shrinking of the excitation wave to ‘zero length’ is a necessary condition, and therefore a ‘cause’ of the block of propagation of excitation waves (Weiss et al. 2000). Comparison of figure 7a,b shows that this is far from true for ionic cardiac models, where such reduction to zero length happens, but only as a very distant consequence, rather than a cause, of the propagation block. The true reason for the block is the disappearance of the fast Na current at the front, observed phenomenologically as its dissipation.

We expect that the condition Embedded Image can also serve as a condition of propagation in the non-stationary situation on the slow space/time scale. Moreover, we conjecture that the dissipation of the front will happen where and when the front runs into a region with Embedded Image. This is illustrated by a simulation shown in figure 7c. It is a solution of the caricature system (4.2), where the excitability parameter J has been maintained slightly above the threshold Embedded Image outside the block domain, and slightly below it within the block domain. As a result, the front propagation has been stopped and never resumed even after the block has been removed. A similar simulation for the quantitatively more accurate fast subsystem of the CRN model, (4.1), is shown in figure 7d. Both agree with what happens in the full model in figure 7a, and both confirm that the condition Embedded Image is relevant for causing front dissipation.

(c) Breakups and self-terminations of re-entrant waves

In two spatial dimensions, the condition Embedded Image may happen to a piece of a wavefront rather than the whole of it. Then instead of a complete block we observe a local block and breakup of the excitation wave. This happened in the episode shown in figure 8.

Figure 8

Analysis of a breakup of a re-entrant wave in a simulation similar to figure 1. Top row: snapshots of the distribution of the transmembrane voltage, at the selected moments of time (designated above the panels). The other three rows: profiles of the key dynamic variables (designated on the left) along the dotted line shown on the top row panels, at the same moments of time. The scale of E is Embedded Image. The scale of Embedded Image is Embedded Image. The scale of j is Embedded Image. Horizontal dash-dotted line on j panels represents Embedded Image. See also Biktashev & Biktasheva (2005).

The white dotted horizontal line on the top panels goes across the region where the propagating wave has been blocked and front has dissipated. The details of how it happened are analysed on the lower three rows, showing profiles of relevant variables along this dotted white line. The second row shows the profiles of E, which lose the sharp front gradient after Embedded Image. The third row shows the peaks of the spatial distribution of the product Embedded Image; the sharpness of these peaks corresponds to the sharp localization of Embedded Image at the front, and their decay accompanies the process of the front dissipation. The most instructive is evolution of the profile of the j variable shown on the bottom row. Consider the column Embedded Image. The gradient of j ahead of the front, i.e. to the left of the peak of Embedded Image, is positive, and the front is moving leftwards, i.e. towards smaller values of j. That is, the front moves into a less excitable area, left there after the previous rotation of the spiral wave. To the right of the peak of Embedded Image the gradient of j is negative which corresponds to the fact that j decreases during the plateau of the AP. Thus, its maximal value at this particular time is observed at the front. This maximal value is, therefore, the value that should be considered in the condition of the dissipation, Embedded Image.

As soon as the front has dissipated (Embedded Image), the profile of j starts to rise, so the maximum of the j profile observed at Embedded Image is the lowest one. From the fact that dissipation has started, we conclude that this maximum is below the critical value Embedded Image. Assuming that dissipation usually happens soon after the condition Embedded Image is satisfied (simulations of (4.1) show that this happens within a few milliseconds), this smallest maximum value should be close to Embedded Image, which gives an empirical method of determining Embedded Image from numerical simulations of complete PDE models. For this particular episode, the empirical value of Embedded Image was found to be approximately 0.3. This is about 50% higher than the Embedded Image predicted for the same range of voltages by (4.1); we attribute this to the approximation Embedded Image, which caused similar errors in the upstroke height and front propagation speeds.

5. Conclusion

Our new asymptotic approach for cardiac excitability equations has significant advantages over the traditional approaches. The fast subsystem, represented by equations (3.3) and (4.1), appears to be typical for cardiac models. This predicts that front propagation cannot happen at a speed slower than a certain minimum and at an excitability parameter lower than a certain minimum. When these conditions are violated the front dissipates and does not recover even after excitability is restored. We have obtained a condition for front dissipation in terms of an inequality involving pre-front values of j and E. This condition can be used for the analysis of breakup and self-termination of re-entrant waves in two and three spatial dimensions.

Editors' note

Please see also related communications in this focussed issue by Benson et al. (2006) and Steinberg et al. (2006).


This work was supported in part by grants from EPSRC (GR/S43498/01, GR/S75314/01) and by an RDF grant from University of Liverpool.


  • One contribution of 13 to a Theme Issue ‘Biomathematical modelling I’.


View Abstract