## Abstract

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

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,(2.1)where *ϵ*> 0 is small, *t* is the ‘slow time’, *T* is the ‘fast time’, , ** X** is the vector of ‘slow variables’ and

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

*Y*We consider Hodgkin & Huxley (1952) (henceforth referred to as HH) and Noble (1962) (henceforth N62) models. Both can be written in the same form as(2.2)where is the total transmembrane current (μA cm^{−2}), *t* is time (ms), *E* is the transmembrane voltage (mV), , are the reversal potentials of sodium, potassium and leakage currents, respectively (mV), are the corresponding maximal specific conductances (mS cm^{−2}), *n*, *m*, *h* are dimensionless ‘gating’ variables, is the specific membrane capacitance (μF cm^{−2}), are the gates' instantaneous equilibrium ‘quasi-stationary’ values and 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, . For a system of differential equations , , we define . The *τ*'s obtained for *m*, *h* and *n* in this way coincide with 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 2*a* 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 remains faster than the pair . This suggests the introduction into system (2.2) of a parameter *ϵ*, which in the limit makes variables *E* and *m* much faster than *n* and *h*:(2.3)

The properties of this system in the limit are shown in figure 2*b*,*c*. The fast transient, corresponding to the AP upstroke, can be studied by changing the independent time variable in (2.3) to and then considering the limit , which gives the *fast subsystem* of two equations for *m* and *E*,(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 2*b*. 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 at all values of *n* and *h* form a two-dimensional *slow manifold* in the four-dimensional phase space of (2.3) with coordinates . Projection of this two-dimensional manifold into the three-dimensional space with coordinates is depicted in figure 2*c*. 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 or equivalently . On the time-scale , 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 can be described by the limit in (2.3), which gives a system of two finite equations and two differential equations,(2.5)The finite equations define the slow manifold and the differential equations define the movement along it.

Figure 2*c* 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 . 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 2*d*. 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 3*a* 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, to describe the difference between the fast and superfast time-scales, and to describe the difference between the slow and fast time-scales. System (2.2) then takes the form(2.6)Consider first the limit . 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, is always close to its quasi-stationary value . Thus, with an error of the order of , we may approximate *m* by and discard the first equation, i.e. *adiabatically eliminate* superfast variable *m*.

With the remaining system of three differential equations for *E*, *h* and *n*, we consider the change of the time variable as before, , and proceed to the limit . This produces the fast subsystem in the form(2.7)Figure 3*b* 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 . Its projection on the plane is shown in figure 3*c*. 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 . 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 , 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 4*a*. 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 in front of it.

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 , is always much smaller than 1. This happens because and is very small for *E* below a certain threshold voltage , and is very small for *E* above another threshold voltage and . So the ranges of almost complete closure of and 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, for , for and , 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 and for the corresponding small parameters. The small parameters should also ensure that and are small in some ranges of *E* but not in others. We identify this smallness with rather than , as it is to compensate the large value of outside the upstroke, and the large value of is described by . We denote the -dependent versions of and as and . In this way, we obtain(3.1)where is the Heaviside function.

The limit gives adiabatic elimination of the *m* gate, .

The analysis of the limit is complicated by a feature incidental to N62 and not found in other cardiac excitation models. The small conductivity of the window current, , is multiplied by a large factor . In N62, the resulting window component of 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 -dependent part of (3.1) can, in the limit , be replaced with a ‘modified N62’ model:(3.2)where , and as before for , for and .

The limit of (3.2) in the fast time gives the fast subsystem(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 4*b*. The horizontal isocline (the red set) is not just a curve but contains a whole domain . The vertical isocline (the blue set) lies entirely within the red set, so the whole line consists of equilibria. An upstroke trajectory may end up in any of the equilibria above , 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 depends on details of approximating function , 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 . The collection of such lines makes a two-dimensional manifold in the three-dimensional space with coordinates . So the fast variable *h* can be adiabatically eliminated on the time-scale . Thus, the slow subsystem, i.e. the limit in (3.2), is(3.4)The phase portrait of this system is shown in figure 4*c*. 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 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 5*a*). 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*, and *w*. As before, we denote the associated small parameter .

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*, and *d* as fast, with associated small parameter .

Similar to N62, the transmembrane voltage is -fast only during the AP upstroke due to -large values of during that period, and is slow otherwise. This is due to nearly perfect switch behaviour of the gates *m* and *h*. The definition of 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:(3.5)where is the fast Na current and represents all other currents. Here we have shown only equations that contain or . All other equations are the same as in the original model.

As before, we adiabatically eliminate the superfast variables in the limit , and turn in the fast time-scale . This gives the fast subsystem(3.6)Note that in (3.6) the equations for and *d* depend on *E*, but equations for *E* and *h* do not depend either on 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, . The phase portrait of (3.6), figure 5*b*, is similar to that of N62, figure 4*b*. 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 5*c* for the shape of the upstroke in the fast time and figure 5*d* 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 and , except for the upstroke: e.g. the peak voltage is overestimated by about 13 mV. This is mostly due to the limit , i.e. replacement of *m* with . So the accuracy could be significantly and easily improved by retracting the limit , which amounts to inclusion of the evolution equation for *m* instead of the finite equation . Most of the qualitative analysis remains valid. However, here for simplicity, we stick to the less accurate but simpler case . Note that in figure 5*b–d*, we took , ; the error introduced by that was small compared to other errors, particularly the error introduced by .

## 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*:(4.1)where and we have put . We consider solutions in the form of propagating fronts. For definiteness, let us assume the fronts propagating leftwards, so , , , , , and . In this formulation, we have three constants characterizing the solutions, the pre-front voltage , the post-front voltage 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 and in it with constants(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 6*a* illustrates a typical front solution. It exists if speed *c* and pre-front voltage satisfy a finite equation involving also the constants *J* and . The resulting dependence of the conduction velocity *c* on excitability *J* for a few selected values of is shown in figure 6*b*. These front solutions exist only for *J* at or above a certain minimum which depends on . For , 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).

The replacement of functions and 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 6*c*. Not only is the topology of the solution set the same, but the overall behaviour of in (4.1) is quite similar to that of 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*, and *w* and the fast variables *h*, 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 , there is a certain minimum excitability and corresponding minimum propagation speed , and for , 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 6*d*: 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 §4*b*, 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 §4*a*, 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 , and the slow inactivation gate *j* of the fast Na current. This dependence gives an equation of motion for the front coordinate ,(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 remains defined, i.e. while . If the front runs into a region where this is not satisfied, its propagation becomes unsustainable.

What will happen then is illustrated in figure 7*a*, where the parameter was varied in space and time. To make the effect more prominent, we did not use smooth variation, but put 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.

This is completely different from the behaviour of a FHN system in similar circumstances, shown in figure 7*b*. 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 7*a*,*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 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 . This is illustrated by a simulation shown in figure 7*c*. It is a solution of the caricature system (4.2), where the excitability parameter *J* has been maintained slightly above the threshold 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 7*d*. Both agree with what happens in the full model in figure 7*a*, and both confirm that the condition is relevant for causing front dissipation.

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

In two spatial dimensions, the condition 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.

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 . The third row shows the peaks of the spatial distribution of the product ; the sharpness of these peaks corresponds to the sharp localization of 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 . The gradient of *j* ahead of the front, i.e. to the left of the peak of , 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 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, .

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

## Acknowledgments

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.

## Footnotes

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

- © 2006 The Royal Society