Originally published as Biophys J. BioFAST on January 13, 2006.
doi:10.1529/biophysj.105.072637
Biophysical Journal 90:2258-2269 (2006)
© 2006 The Biophysical Society
Conditions for Propagation and Block of Excitation in an Asymptotic Model of Atrial Tissue
Radostin D. Simitev and
Vadim N. Biktashev
Department of Mathematical Sciences, University of Liverpool, Liverpool, United Kingdom
Correspondence: Address reprint requests to Vadim N. Biktashev, Dept. of Mathematical Sciences, University of Liverpool, Liverpool L69 7ZL, UK. Tel.: 44-151-7944004; Fax: 44-151-7944061; E-mail: vnb{at}liv.ac.uk.
 |
ABSTRACT
|
|---|
Detailed ionic models of cardiac cells are difficult for numerical simulations because they consist of a large number of equations and contain small parameters. The presence of small parameters, however, may be used for asymptotic reduction of the models. Earlier results have shown that the asymptotics of cardiac equations are nonstandard. Here we apply such a novel asymptotic method to an ionic model of human atrial tissue to obtain a reduced but accurate model for the description of excitation fronts. Numerical simulations of spiral waves in atrial tissue show that wave fronts of propagating action potentials break up and self-terminate. Our model, in particular, yields a simple analytical criterion of propagation block, which is similar in purpose but completely different in nature to the "Maxwell rule" in the FitzHugh-Nagumo type models. Our new criterion agrees with direct numerical simulations of breakup of reentrant waves.
|
INTRODUCTION
|
|---|
Refractoriness is a fundamental characteristic of biological excitable media, including cardiac tissues. The boundary between absolute and relative refractoriness can be defined as the boundary between the ability and the inability of the medium to conduct excitation waves (1
). Transient conduction block is thought to be a key event in the initiation of reentrant arrhythmias and in the development and the self-perpetuation of atrial and ventricular fibrillation (2–5). So it is important to understand well the immediate causes and conditions of propagation blocks and sudden breakups in such nonstationary regimes. The aim of this work is to improve this understanding via analysis of a mathematical model of human atrial tissue (6).
Kohl et al. (7) distinguish two types of single-cell cardiac models: "membrane potential models" and "ionic current models". The membrane potential models attempt to represent cellular electrical activity by describing, with a minimal number of equations, the spatio-temporal course of changes in membrane potential. Their equations are constructed using dynamical systems arguments to caricature various properties and processes of cardiac function. Examples of this type of models start with the mathematical description of heartbeat as a relaxation oscillator by van der Pol and van der Mark (8) and continue to play an important role in describing biophysical behavior (9) with the most successful one arguably being the FitzHugh-Nagumo equations (10,11),
 | (1) |
where V and g are dynamical variables corresponding to the action potential and the cardiac current gating variables, 
V, g, 
, and ß are parameters, and D is a diffusion constant. Further examples of such models can be found in Aliev and Panfilov, Pertsov and Panfilov, Barkley, and Winfree (12–15), among others. An attractive feature of this approach is that, along with a reasonable description of excitability, threshold, plateau, and refractoriness, it focuses on generic equations that can often be treated analytically and their dynamical properties can be extended and applied to very different physical, chemical, or biological problems of similar mathematical structure. The main drawback of these models, however, is their lack of an explicit correspondence between model components and constituent parts of the biological system, e.g., ion channels and transporter proteins. The second type of models, the ionic current models, attempt to model action potential (AP) behavior on the basis of ion fluxes in as much detail as possible to fit experimental data and predict behavior under previously untested conditions. A major breakthrough in this direction of cell modeling was the work of Hodgkin and Huxley (16), representing the first complete quantitative description of the giant squid axon. The ionic concept was applied to cardiac cells by Noble (17,18) and there are now ionic models of sinoatrial node pacemaker cells, e.g. (19); atrial myocytes, e.g. (20); Purkinje fibers, e.g. (21); ventricular myocytes, e.g. (22,23); and cardiac connective tissue cells, e.g. (24). This is only an incomplete list and the collection of available models continues to expand. The ionic models have been successfully applied to study various conditions of metabolic activity and excitation-contraction coupling, feedback mechanisms, response to drugs, etc. For recent reviews of detailed ionic models, their computational aspects, and applications, we refer to the reviews of Kohl et al. (7) and Clayton (25). However, these models are very complicated and have to be studied mostly numerically. Their numerical study is aggravated by stiffness of the equations, i.e., broad range of characteristic timescales of dynamic variables caused by numerous small parameters of the models.
An attractive compromise is exemplified by the model of Fenton and Karma (26), which combines the simplicity of only three differential equations with realistic description of (crudely) the AP shape and (rather nicely) the dependence of the AP duration and front propagation speed on the diastolic interval, i.e., "restitution curves". Unlike the earlier two-component model by Aliev and Panfilov (12), it has a structure similar to that of true ionic models, and its parameters have been fitted to mimic properties of selected four detailed ventricular myocyte models. It is simpler than later proposed models of the same "intermediate" kind such as Bernus et al. (27). However, this deservedly popular model has not been in any way "derived" from any detailed model, so it is only reliable within the phenomenology on which it has been validated, i.e., normal or premature APs, but not propagation blocks.
The problem of conditions for propagation has an elegant solution for the FitzHugh-Nagumo system Eqs. 1 and its generalizations, within an asymptotic theory exploiting the difference of timescales of different variables, such as g << V in case of Eqs. 1 (28). The answer is formulated in terms of the instantaneous values of the slow variables (g in Eqs. 1), and claims that excitation will propagate if the definite integral of the kinetic term on the right-hand side of the equation for the fast variable (V in Eqs. 1), between the lower and the upper quasi-stationary states, is positive (see Eq. 4.5 in Fife (29)). This is similar to Maxwell's "equal areas" rule in the theory of phase transitions (see section 9.3 in Haken (30)). In case of Eqs. 1, this rule boils down to an inequality for the slow variable g: excitation front will propagate if the value of g at it is less than a certain g*. However, FitzHugh-Nagumo-type models completely misrepresent the idiosyncratic "front dissipation" scenario by which propagation block happens in the ionic current models (31). The reason is that small parameters in such models appear in essentially different ways from the one assumed by the standard asymptotic theory (32,33). So, this elegant "Maxwell rule" solution is not applicable to any realistic models.
We have developed an alternative asymptotic approach based on special mathematical properties of the detailed ionic models, not captured by the standard theory (34). This approach demonstrated excellent quantitative accuracy for APs in isolated Noble-1962 model cells (33), and correctly, on a qualitative level, described the front dissipation mechanism of breakup of reentrant waves in the Courtemanche et al. (6) model of human atrial tissue, although quantitative correspondence with the full model was poor (35). In this article we suggest, for the first time, to our knowledge, a refined simplified asymptotic model of a cardiac excitation front, which provides numerically accurate prediction of the front propagation velocity (within 16%) and its profile (within 0.7 mV). It also gives an analytical condition for propagation block in a reentrant wave, expressed as a simple inequality involving the slow inactivation gate j of the fast sodium current. The condition is in excellent agreement with results of direct numerical simulations of the Courtemanche et al. (6) full ionic model of 21 partial differential equations.
The article is organized as follows. In the next section, we introduce simplified model equations and discuss their properties. Analytical solutions are then presented for a piecewise linear "caricature" version of our simplified model, followed by numerical results and a two-dimensional test. The article concludes with a discussion of results and questions open for future studies.
|
MATHEMATICAL FORMULATION OF THE MODEL EQUATIONS
|
|---|
Asymptotic reduction
In this section, we briefly summarize the asymptotic arguments of Biktasheva et al. (35) relevant to our present purposes. We rewrite the Courtemanche et al. (6) model in the following one-parameter form:
 | (2) |
where D is the voltage diffusion constant, is a small parameter used for the asymptotics, 
is the curvature of the propagating front, 
() is the Heaviside function, 
'I() is the sum of all currents except the fast sodium current INa, the dynamic variables V, m, h, ua, oa, and d are defined in Courtemanche et al. (6), U = (j, oi, ..., Nai, Ki, ...)T is the vector of all other, slower variables, and F is the vector of the corresponding right-hand sides. The rationale for this parameterization is:- The dynamic variables V, m, h, ua, w, oa, and d are "fast variables", i.e., they change significantly during the upstroke of a typical AP potential, unlike all other variables that change only slightly during that period. The relative speed of the dynamical variables is estimated by comparing the magnitude of their corresponding "timescale functions" as shown in Fig. 1 A. For a system of differential equations
the timescale functions are defined as 
i(y) 
|(dFi/dyi)–1|, i = 1...N and coincide with the functions already present in Eqs. 2.
- A specific feature of V is that it is fast only because of one of the terms on the right-hand side, the large current INa, whereas other currents are not that large and so do not have the large coefficient –1 in front of them.
- The fast sodium current INa is only large during the upstroke of the AP, and not that large otherwise as illustrated in Fig. 1 D. This is due to the fact that either gate m or gate h or both are almost closed outside the upstroke since their quasi-stationary values
and 
are small there as seen in Fig. 1 B. Thus in the limit 
0, functions and have to be considered zero in certain overlapping intervals V 
(– 
, Vm] and V [Vh, + ), and Vh 
Vm, hence the representations 
and 
- The term
in the first equation represents the effect of the front curvature for waves propagating in two or three spatial dimensions. Derivation of this term using asymptotic arguments can be found, e.g., in Tyson and Keener (28). A simple rule-of-thumb way to understand it is this. Imagine a circular wave in two spatial dimensions. The diffusion term in the equation for V then has the form 
where R is the polar radius. If R at the front is large, its instant curvature 
changes slowly as the front propagates, and can be replaced with a constant for long time intervals. Considering R as a new X coordinate, we then get Eqs. 2.
These aspects, as applied to the fast sodium current, have been shown to be crucial for the correct description of the propagation block (31). In particular, it is important that the h gate is included among the fast variables. The particular importance of h dynamics at the fringe of excitability has been noted before, e.g., for the modified Beeler-Reuter model (36). A more detailed discussion of the parameterization given by Eqs. 2 can be found in Biktasheva et al. (35).
A change of variables t = –1T, x = (D)–1/2X, 
and subsequently the limit 0 transforms Eqs. 2 into
 | (3) |
In other words, we consider the fast timescale on which the upstroke of the AP happens, neglect the variations of slow variables during this period as well as all transmembrane currents except INa, as they do not make significant contribution during this period and replace 
and 
with zero when they are small. (A change of the value of D is equivalent to rescaling of the spatial coordinate, and is not critical to any of the questions considered here. To operate with dimensional velocity, we assume the value of the diffusion coefficient D = 0.03125 mm2/ms, as in our earlier publications (35,37). Increase of the diffusion coefficient to, say, D = 0.1 mm2/ms, raises the propagation velocity from 0.28 mm/ms in Table 1 to 0.50 mm/ms, in full agreement, e.g., with results of Xie et al. (38) for the same model.)
View this table:
[in this window]
[in a new window]
|
TABLE 1 A comparison of the wave speed C, postfront voltage amplitudes V and the maximum rate of AP rise (dV/ dt)max of various approximations to the Courtemanche et al. (6) model
|
|
In the resulting Eqs. 3, the first three equations for V, m, and h form a closed subsystem. The following four equations for ua, w, oa, and d can be solved if V(x, t) is known but do not affect its dynamics, and the rest of the equations state that all other variables remain unchanged. Hence we concentrate on the first three equations as the system describing propagation of an AP front or its failure. The above derivation procedure does not give a precise definition of the functions H(V) and M(V); it only requires that these are reasonably close to and for those values of V where these functions are not small. Here "reasonably close" means that replacement of with H(V) (Vh – V) and with M(V) (V – Vm) does not change significantly the solutions of interest, i.e., the propagating fronts. We have found that the simplest approximation in the form M(V) = 1, H(V) = 1 works well enough. This is demonstrated in Table 1 where various choices of M(V) and H(V) are tested. So, ultimately, we consider the following system:
 | (4a) |
 | (4b) |
 | (4c) |
where
 | (5a) |
 | (5b) |
All parameters and functions here are defined as in Courtemanche et al. (6) except the new "gate threshold" parameters Vh and Vm, which are chosen from the conditions 
and 
As follows from the derivation, variable j, the slow inactivation gate of the fast sodium current, acts as a parameter of the model. It is the only one of all slow variables included in the vector U that affects our fast subsystem. We say that it describes the "excitability" of the tissue. Notice that it is a multiplier of gNa, so a reduced availability of the fast sodium channels, e.g., as under tetrodotoxin (39) or arguably in Brugada syndrome (40) can be formally described by a reduced value of the parameter j.
Before proceeding to the analysis of the simplified three-variable model defined by Eqs. 4, we wish to demonstrate that it is a good approximation of the full model of Courtemanche et al. (6) both on a qualitative and a quantitative level. On the qualitative level, we show that a temporary obstacle leads to a dissipation of the front. This is illustrated in Fig. 2, which shows propagation of the AP into a region in time and space where the excitability of the tissue is artificially suppressed. The sharp wave fronts of the model of Courtemanche et al. (6) as well as of Eqs. 4 stop propagating and start to spread diffusively once they reach the blocked zone. The propagation does not resume after the block is removed. This behavior is completely different from that of the FitzHugh-Nagumo system of Eqs. 1 in which even though the propagation is blocked for nearly the whole duration of the AP, the wave resumes once the block is removed. Table 1 illustrates, on the quantitative level, the accuracy of Eqs. 4 as an approximation of the full model of Courtemanche et al. (6).
It is a popular concept going back to classical works (e.g., 41) that the fast activation gate m is considered a "fast variable" and is "adiabatically eliminated" since most of the time, except possibly during a very short transient, it is close to its quasi-stationary value 
Hence the model can be simplified by replacing m with and eliminating the equation for m,
 | (6) |
We have explored this possibility for the model of Courtemanche et al. (6) in Biktasheva et al. (35). Equations 6 are qualitatively correct, i.e., they still show front dissipation on collision with a temporary obstacle, but make a large error in the front propagation speed, as demonstrated in Table 1.
Traveling waves and reduction to ordinary differential equation of the three-variable model
To find out when propagation of excitation is possible in our simplified model and when it will be blocked, we study solutions in the form of propagating fronts as well as the conditions of existence of such solutions.
We look for solutions in the form of a front propagating with a constant speed and shape. So we use the ansatz F(z) = F(x + ct) for F = V, h, m, where z = x + ct is a "traveling wave coordinate" and c is the dimensionless wave speed of the front, related to the dimensional speed C by c = (/D)1/2C. Then Eqs. 4 reduce to a system of autonomous ordinary differential equations (ODE),
 | (7a) |
 | (7b) |
 | (7c) |
where the boundary conditions are given by
 | (8a) |
 | (8b) |
 | (8c) |
Here V
and V are the pre- and postfront voltages.
Equations 7 represent a system of fourth order so its general solution depends on four arbitrary constants. Together with constants V, V, and c, this makes seven constants to be determined from the six boundary conditions in Eqs. 8. Thus, we should have a one-parameter family of solutions, i.e., one of the parameters (V, V, c) can be chosen arbitrary from a certain range. A natural choice is V because the prefront voltage acts as an initial condition for a propagating front in the tissue, and because in our study of the conditions for propagation it is most conveniently treated as a parameter rather than as an unknown.
|
ANALYTICAL STUDY OF THE REDUCED MODEL
|
|---|
An exactly solvable caricature model
The parameter-counting arguments given in the previous section make it plausible that the problem defined by Eqs. 7 with boundary conditions of Eqs. 8 has a one-parameter family of traveling wave-front solutions. However, the problem is posed in a highly unusual way since the asymptotic prefront and postfront states are not stable isolated equilibria but belong to continua of equilibria and thus are only neutrally stable. We are not aware of any general theorems that would guarantee existence of solutions of a nonlinear boundary value-eigenvalue problem of this kind. For the two-component model of Eqs. 6 considered in Biktasheva et al. (35), this worry has been alleviated by the fact that there is a "caricature" model, which has the same structure as Eqs. 6, including the structure and stability of the equilibrium set and which admits an exact and exhaustive analytical study (31). Fortunately, a similar "caricature" exists for our present three-variable problem as well. We replace functions 
h(V), and m(V) defined in Eqs. 5 with constants. The choice of the constants is somewhat arbitrary. We assume that the events in the beginning of the interval z [
, + ), where V is just above Vm, are most important for the front propagation. So for numerical illustrations we choose the values of constants 
h, and m as the values of the corresponding functions in Eqs. 5 at some fixed value of the voltage V. We set the z axis so that V(0) = Vh, and then V() = Vm for some > 0 still to be determined. We demand that the solutions for the unknowns V, h, and m are continuous and that V is smooth at the internal boundary points.
In this formulation, Eqs. 7b and 7c decouple from Eq. 7a and from each other and are solved separately. The solutions of these first-order linear ODE with constant coefficients are given by Eqs. 10b and 10c, respectively. It follows that in the interval V Vm, Eq. 7 is a linear homogeneous ODE with constant coefficients, and its solution given at the first row of Eq. 10a satisfies the boundary conditions V(– ) = V, V(0) = Vh, and V() = Vm, provided that the internal boundary point is given by Eq. 12. To solve the linear inhomogeneous Eq. 10 in the interval V 
Vm, we note that its inhomogeneous term 
is a sum of exponentials
 | (9) |
and terms proportional to nh will appear in the solution due to the expression for Bn. Imposing the boundary conditions at the internal point V() = Vm and at infinity V() = V, we obtain the solution in this interval given at the second row of Eq. 10a. Finally, the wave speed c is fixed by Eq. 11b from the requirement that the solution for V(z) is smooth at the internal boundary point . To summarize, the solution of Eqs. 7 and 8 is
 | (10a) |
 | (10b) |
 | (10c) |
where the prefront voltage V, the postfront voltage V, and the wave speed c are related by
 | (11a) |
 | (11b) |
the distance between points V = Vh and V = Vm is
 | (12) |
and An(c, z) and an(c) are abbreviations for
 | (13a) |
 | (13b) |
In the limit m 0, this solution tends to the solution of the two-component model of (42), as expected.
The accurate expression in Eq. 5a for the sodium current 
vanishes for V = VNa, which, in particular, means that the transmembrane voltage never exceeds VNa. So, replacing this function with a constant changes the properties of the system qualitatively. Even bigger discrepancies are expected to occur from replacing the h(V) and m(V) by constants because these functions vary by an order of magnitude in the range between the pre- and the postfront voltage. It is surprising, however, that even this rough approximation produces results that, with exception of the postfront voltage, are within several percent from the solution of the detailed ionic model (6) and certainly capture its qualitative features as can be seen in Fig. 3, where the constants are chosen at V = Vm, i.e., 
h(Vm), and m(Vm). This relatively good agreement is not due to this special choice of parameter values. Indeed, the caricature model and its solution Eqs. 10 involve the parameters h, m, 
, V, and j. The dependence on the curvature is negligible in comparison to the deviation of the solution Eqs. 10 of the caricature model from the numerical solution of the three-variable model Eqs. 10. The dependence on the prefront voltage V and the excitability parameter j is discussed in the subsection immediately following and represented in Figs. 4 and 6. The parameters h, and m, on the other hand, are somewhat arbitrary but to achieve a good agreement with the original system given by Eqs. 7, we choose these values as the values of the corresponding functions in Eqs. 5 at various values of V. In Fig. 4, the relationship between the wave speed c and the excitation parameter j for several such choices of V is presented. It can be seen that such a variation of the values of h, and m does not lead to significant qualitative changes in the solution Eqs. 10 of the caricature model. Figs. 3 and 4 also show, for comparison, the numerical solutions of the detailed ionic model of Courtemanche et al. (6) and of the full three-variable model of Eqs. 7, which will be described in detail in the next section.
The condition for propagation
Equation 11b defines c as a smooth function of the parameters within a certain domain. The boundary of this domain is associated with the propagation failure. Not all parameters, h, m, , V, and j, entering Eq. 11b are of equal importance. We consider here = 0 and postpone the investigation of the effects of curvature to the next section. Parameters h, and m represent well-defined properties of the tissue, albeit changeable depending on physiological conditions. On the other hand, parameters j and V are not model constants, but "slowly varying" dynamic quantities: j remains approximately constant throughout the front, and V represents the transmembrane voltage ahead of the front, but both can vary widely on large scales between different fronts. Hence we need to determine the singular points of the dispersion relation in Eq. 11b with respect to j and V. Similarly to the two-component caricature (31), Eq. 11b is a transcendental equation for c, but it is easily solvable for the excitation parameter j:
 | (14) |
The resulting relationship of j and c for a selected value of V is shown in Fig. 4. This figure reveals a bifurcation. For values of j lower than some jmin, no traveling wave solutions exist. After a bifurcation at j > jmin, two solutions with different speeds are possible. Our direct numerical simulations of Eqs. 4 as well as studies of the two-component caricature model by Hinch (43) suggest that the solutions of the lower branch are unstable. The bifurcation point jmin can be determined from the condition that j has a minimum with respect to c at this point and therefore satisfies
 | (15) |
This produces, with j(c) defined by Eq. 14, a quintic polynomial equation for c2.
Activation of the sodium current is possible because m = h, permitting transient channel opening and current flow through the cell membrane. The ratio h/m is a function of V in the full model, and is a constant in Eqs. 7. The minimal value of this ratio, necessary for propagation, is shown in Fig. 5 as a function of various choices of m, and j; it is obtained by numerical solution of the algebraic equation Eq. 11b. The smallness of m/h allows approximate solution of the above mentioned quintic equation for c2. We set
 | (16) |
Substituting this expansion in Eq. 15 and discarding the small terms of order O(m) gives the zeroth-order approximation to the solution as a function of the prefront voltage V:
 | (17) |
This limit corresponds to the two-variable caricature (31). For any given value of the prefront voltage, the value of j must be larger than jmin for wave fronts to propagate. Although lacking sufficient accuracy, the zeroth-order approximation given by Eq. 17 reproduces qualitatively well the conditions for propagation and dissipation of excitation fronts in the model of Courtemanche et al. (6). Analogously, discarding the small terms of order 
gives the first-order approximation,
 | (18) |
This approximation is already very good and changes insignificantly as more terms are considered in Eq. 16 (see Fig. 6).
|
NUMERICAL RESULTS
|
|---|
Propagating front solutions
We solved Eqs. 7–8 numerically, using the method described in the Appendix. The results are shown in Figs. 3, 4, and 7. Fig. 3 offers a comparison of the shapes of the solution of Eqs. 7 with a snapshot of a traveling wave solution of the full model of Courtemanche et al. (6). The values of the wave speed and the postfront voltage are presented in Table 1 and also show an excellent agreement. This confirms our assumptions that the fronts of traveling waves in the full model have constant speed and shape and thus satisfy an ODE system, and that j remains approximately constant during the front. Fig. 7 shows the wave speed c as a function of two of the parameters of the problem, the prefront voltage V and the excitability parameter j. For every value of j and V from a certain domain, two values of the wave speed c are possible, which is similar to the solutions of the caricature model. The smaller values of c are not observed in the partial differential equations simulation of the full model. This is a strong indication that they are unstable.
View larger version (27K):
[in this window]
[in a new window]
|
FIGURE 7 Wave speed C as a function of j and V, for the model of Eqs. 7. Rapid changes are indicated by a higher density of curves. The thick dotted line on the base represents the threshold value jmin and may be compared to the results in Fig. 6.
|
|
The condition for propagation
In this subsection, we report numerical values for the threshold of excitability jmin below which wave fronts are not sustainable and have to dissipate, as predicted by the reduced three-variable model of Eqs. 7–8. Fig. 8 presents jmin as a function of the prefront voltage V. The curve jmin(V) represents a boundary in the space of the slow variables (V, j) which separates the region of relative refractoriness where excitation fronts are possible, even though possibly slowed down, from the region of absolute refractoriness where excitation fronts cannot propagate at all. In practice, however, we can reduce the condition of the absolute refractoriness even further. This is possible because typical APs have their tails very closely following one path on the (V, j) plane. This property is known for cardiac models; e.g., Vinet and Roberge (36) present an evidence for the modified Beeler-Reuter model that the dynamics of recovery from an AP do not depend on details of how that AP has been initiated. Therefore of the whole curve (V, jmin(V)), only one point is important—its intersection with the curve (V(t), j(t)), representing a typical AP tail. For the model Courtemanche et al. (6) considered here, we simply state the existence of this universal (V(t), j(t)) curve as an "experimental fact". This is illustrated in Fig. 8, where we plot the curve (V, jmin(V)) together with projections of a selected set of AP trajectories. The AP solutions were obtained for a space-clamped version of Courtemanche et al. (6) with initial conditions for j and V as shown in the figure and all other variables in their resting states. These trajectories allow us to follow the correlation between the transient of j and the AP V. Indeed, in the tail of an AP solution, the curve j versus V is almost independent of the way the AP is initiated. As a result, the projections of the trajectories (V(t), j(t)) intersect the critical curve (V, jmin(V)) in a small vicinity of one point, (j*, V*) = (0.2975 ± 0.0015, –72.5 ± 0.5). This result suggests the following interpretation. As a wave front propagating into the tail of a preceding wave reaches a point in the state corresponding to this "absolute refractoriness" point (j*, V*), it will stop because of insufficient excitability of the medium, and dissipate.
In a broader context, in the front propagation speed, c is a function of j and V in the relative refractoriness region of the (V, j) plane, so the highly correlated dependencies of V(t) and j(t) in the wake of an AP mean that c at a particular point becomes a fixed function of time. This makes it possible to describe c in terms of the diastolic interval DI, i.e., the time passed after the end of the preceding AP. This dependence, known as dispersion curve or velocity restitution curve, is an important tool in simplified analysis of complex regimes of excitation propagation (44–48).
Propagation block in two dimensions
In two spatial dimensions, the condition of dissipation j < j* may happen to a piece of a wave front rather than the whole of it. In that case, we observe a local block and a breakup of the excitation wave. Fig. 9 shows how it happens in a two-dimensional simulation of the detailed model of Courtemanche et al. (6). A spiral wave was initiated by a cross-field protocol. This spiral wave develops instability, breaks up from time to time, and eventually self-terminates. This is one of the simulations discussed in detail in Biktasheva et al. (35). Here we use it to test our newly obtained criterion of propagation block. The red color component represents the V field, white for the resting state, and maximum for the AP peak. This is superimposed onto an all-or-none representation of the j field, with white for j > j* and blue for j j*. Thus the red rim represents the "active front" zone where excitation has already happened but j gates are not deactivated yet; most of the excited region is in shades of purple representing the gradual decay of the AP with j deactivated. The wave ends up with a blue tail, which corresponds to V already close to the resting potential but j not yet recovered and still below j*. So the blue zone is where there is no excitation, but propagation of excitation wave is impossible, i.e., absolutely refractory zone. The white zone after the tail and before the new front is therefore relative refractory zone, where front propagation is possible. Thus, in terms of the color coding of Fig. 9, the prediction of the theory is: the wave front will be blocked and dissipate where and when it reaches the blue zone, and only there and then. This is exactly what happens in the shown panels: the red front touches the blue tail, first at the third panel, at the point indicated by the white arrow, and subsequently in its vicinity. The excitation front stops in that vicinity and dissipates. So we have a breakup of the front.
View larger version (77K):
[in this window]
[in a new window]
|
FIGURE 9 Local propagation block, dissipation, and breakup of the front of a reentrant excitation wave. The density plots represent the distribution of the transmembrane voltage V (red component) in regions of superthreshold (white) and of subthreshold (blue) excitability j. The white arrow indicates the time and place the propagation block begins. The time increases from A to F with  t = 20 ms; size of the simulation domain is 75 mm x 75 mm.
|
|
The analysis of the numerics, which ran for the total of 7400 ms until self-termination of the spiral and showed four episodes of front breakup, has confirmed that in all cases the breakup happened if and only if the front reached the blue region j j*.
Curvature effects
Since we attempt to compare the results of our one-dimensional model to simulations of spiral waves in two-dimensions, it is important to explore the dependence of the solution on the curvature of the front. The standard theory says that in two dimensions, the normal velocity of the wave front needs to be corrected by the term 
where is the typical width of the wave front (28). The speed-curvature diagram presented in Fig. 10 A shows that in our simplified model, this relationship is satisfied to rather large values of | |. Our choice of boundary conditions in Eqs. 8 assumes that the excitation fronts propagate from right to left, so positive values of the curvature correspond to concave fronts. Only at very small values of the radius of curvature of the order of 0.3 mm for j = 1 the wave speed shows a nonlinear dependence on curvature as seen in the inset to Fig. 10 B. This part of the figure also demonstrates that there is a critical value of the curvature for which the excitation wave stops to propagate as well as an unstable branch of the solution. However, these phenomena occur at very large curvatures that are far outside of the range of values of | | < 0.1 mm–1 observed in the two-dimensional simulations of Fig. 9.
The most important question with respect to our study is whether the curvature changes significantly the critical value of the excitation parameter j* below which the wave fronts fail to propagate. To answer this question, we present Fig. 10 C in which the wave speed c is shown as a function of j for three values of the curvature corresponding to a noncurved front and to convex and concave fronts with radius of curvature equal to 10 mm. The values of jmin for these three cases differ only slightly. So, the propagation blocks in our simulations do not depend significantly on the curvature of the front.
This conclusion is valid for the particular cardiac model (6) and for the particular context. In Comtois and Vinet (49), the minimal diastolic interval, defined as time from the moment V = –50 mV to the moment propagation becomes possible again, depended only slightly on curvature for the modified Beeler-Reuter model at standard parameters, but was much more pronounced when j was artificially increased sixfold. The simplest explanation of this difference is that the small variation of jmin due to the curvature takes much longer for j(t) to make if 
j/t is very small, so even that small variation jmin becomes significant.
|
CONCLUSIONS
|
|---|
In this article, we have shown that propagation of excitation and its block in the Courtemanche et al. (6) model of human atrial tissue can be successfully predicted by a simplified model of the excitation front, obtained by an asymptotic description focused on the fast sodium current, 
Whereas it was known that main qualitative features of the INa-driven fronts can be described by a two-component model for V and h, we have now found that for good quantitative predictions, one must also take into account the dynamics of m gates. Thus, we have proposed a three-component description of the propagating excitation fronts given by Eqs. 4. We have obtained an exact analytical solution for a piecewise-linear "caricature" three-component model of Eqs. 4. For an appropriate choice of parameters, it reproduces the key qualitative features of the accurate three-component model of Eqs. 4 and gives a correct order of magnitude quantitatively. Numerical solution of the automodel equation of the proposed three-component model of Eqs. 4 gives a very accurate prediction of propagation block in two-dimensional reentrant waves. For the given model, this reduces to a condition involving the prefront values of V and j, or even in terms of j alone. This provides the sought-for operational definition of absolute refractoriness in terms of j, simple and efficient.
The success of the propagation block prediction justifies the assumptions made on the asymptotic structure, i.e., appearance of the small parameter of Eqs. 2, and also confirms that two-dimensional effects, e.g., front curvature, do not significantly affect the propagation block conditions, at least in the particular simulation.
As the description and role of INa are fairly universal in cardiac models, most of the results should be applicable to other models. However, some other cardiac models may require a more complicated description. For instance, the contemporary "Markovian" description of INa (e.g., (50)) is very different from the classical m3jh scheme. Also, propagation in ventricular tissue in certain circumstances can be essentially supported by L-type calcium current rather than mostly INa alone (51).
|
APPENDIX: NUMERICAL METHOD
|
|---|
For a numerical solution, the problem needs to be formulated on a finite interval z [zmin, zmax] rather than on the open interval z (– , ). Furthermore, because of the piecewise definition of the problem, this interval must be separated in three parts—[zmin, 0], [0, ], and [, zmax] as discussed in section "Analytical study of the reduced model". The standard numerical methods we use require that the problem is posed on a single interval, for instance y [0, L]. So we use the mapping
 | (19) |
to transform Eqs. 7 as follows:
 | (20) |
where the subscripts 1, 2, and 3 denote the variables corresponding to the three subintervals. Here, the end of the second subinterval is an unknown parameter and together with the wave speed c and the postfront voltage V must be determined as a part of the solution. Because these unknowns are constants, their derivatives must vanish, which leads to the introduction of the last three equations in Eqs. 20.
The boundary conditions in Eqs. 8 at infinity are substituted by
 | (21) |
where u is the vector of unknown variables and v is a vector of small perturbations, obtained as a solution of Eqs. 7 linearized about Eqs. 8. Together with the implicit assumptions V(0) = Vm and V() = Vh, which break the translational invariance and the additional requirements that the solutions must be continuous functions of z and that V(z) must be smooth, the necessary 15 conditions are
 | (22) |
We use the boundary-value problem solver D02RAF of the Numerical Algorithms Group numerical library, which employs a finite-difference discretization coupled to a deferred correction technique and Newton iteration (52). The analytical solution given in Eqs. 10 is used as an initial approximation to start the correction process. The method proves to be very robust over a large range of parameters.
|
ACKNOWLEDGEMENTS
|
|---|
The authors are grateful to I. V. Biktasheva for sharing her experience of simulation of the Courtermanche et al. model (6), to H. Zhang and P. Hunter for inspiring discussions related to this manuscript, and to the anonymous referees for constructive criticism and helpful suggestions.
This work is supported by Engineering and Physical Sciences Research Council grants GR/S43498/01 and GR/S75314/01.
Submitted on August 14, 2005;
accepted for publication December 12, 2005.
|
REFERENCES
|
|---|
1. Krinsky, V. I. 1966. Spread of excitation in an inhomogeneous medium (state similar to cardiac fibrillation). Biofizika. 11:776–784.2. Moe, G. K. 1962. On the multiple wavelet hypothesis of atrial fibrillation. Arch. Int. Pharmacodyn. Ther. 140:183–188.
3. Weiss, J. N., P. S. Chen, Z. Qu, H. S. Karagueuzian, and A. Garfinkel. 2000. Ventricular fibrillation: How do we stop the waves from breaking? Circ. Res. 87:1103–1107.[Abstract/Free Full Text]
4. Panfilov, A., and A. Pertsov. 2001. Ventricular fibrillation: evolution of the multiple-wavelet hypothesis. Philos. Trans. R. Soc. Lond. A. 359:1315–1325.[CrossRef]
5. Kléber, A. G., and Y. Rudy. 2004. Basic mechanisms of cardiac impulse propagation and associated
TOP
ABSTRACT
INTRODUCTION
(C) Transmembrane voltage V as a function of time. (D) Main ionic currents versus time. Iin = Ib,Na + INaK + ICa,L + Ib,Ca + INaCa and Iout = Ip,Ca + IK1 + Ito + IKur + IKr + IKs + Ib,K are the sums of all inward and outward currents, respectively, and the individual currents are described in Courtemanche et al. (6). The results are obtained for a space-clamped version of the model at values of the parameters as given in Courtemanche et al. (6). In C and D, a typical AP is triggered by initializing the transmembrane voltage to a nonequilibrium value of V = –20 mV.
) in the models of (A) Courtemanche et al. (6
The solution of the model of Courtemanche et al. (6
is indicated by a dash-dotted line.
mm–1.
are fixed to the values of the corresponding functions in Eqs. 5 at a selected voltage V = V0, the prefront voltage is V