Models of unidirectional propagation in heterogeneous excitable media

Models of unidirectional propagation in heterogeneous excitable media

Applied Mathematics and Computation 216 (2010) 1337–1348 Contents lists available at ScienceDirect Applied Mathematics and Computation journal homep...

540KB Sizes 2 Downloads 83 Views

Applied Mathematics and Computation 216 (2010) 1337–1348

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage:

Models of unidirectional propagation in heterogeneous excitable media John G. Alford Department of Mathematics and Statistics, Sam Houston State University, P.O. Box 2206, Huntsville, TX 77341-2206, United States

a r t i c l e

i n f o

Keywords: Differential equation models Unidirectional propagation Excitable media Heterogeneous parameters

a b s t r a c t Two differential equation models of excitable media (threshold and recovery kinetics) with solutions that exhibit unidirectional propagation are presented. It is shown that unidirectional propagation in heterogeneous excitable media with non-oscillatory kinetics can be initiated from homogeneous initial data. Simulations on a reaction–diffusion model with FitzHugh–Nagumo kinetics and spatially heterogeneous parameters yields a rotating wave on a one-dimensional circular spatial domain. An ordinary differential equation model with four semi-coupled excitable cells and heterogeneous parameters is analyzed to determine a critical parameter region over which unidirectional propagation may occur. Ó 2010 Elsevier Inc. All rights reserved.

1. Introduction The nerve and cardiac cells are examples of excitable cells. If a stimulus of sufficient magnitude is applied to an excitable cell for a brief time interval, the cell membrane voltage will reach a threshold value and exhibit a large deviation away from its equilibrium voltage before returning to rest. This is known as an action potential. If consecutive stimuli are applied, an excitable cell requires time between stimuli to recover and initiate another action potential. The dynamics of threshold and recovery are fundamental to the mathematical models presented here. An action potential may propagate from one cell to another when excitable cells are coupled or spatially distributed. Propagation of the action potential is effected by heterogeneities in the medium such as variations in intercellular coupling strength, thresholds of excitability, resting potential, and stimulus strength. When the cell-to-cell spread of an action potential is blocked, there is propagation failure. Propagation failure in spatially discrete and continuous, inhomogeneous excitable media has been widely investigated, e.g. see [1–8]. A wave of excitation which circulates in the heart around a closed loop is called a reentrant arrhythmia [9]. The onedimensional (ring) model of such an arrhythmia is considered in [10–13]. Unidirectional block must occur in order to establish a reentrant arrhythmia on a ring since bidirectional propagation leads to annihilation of oppositely traveling wavefronts. In [14] we investigated the bifurcation of non-stationary steady-states on a ring (rotating waves) in excitable media without heterogeneities using the FitzHugh–Nagumo (FHN) [15,16] reaction–diffusion model [10]. Simulations of the bifurcating solutions used a rotating wave profile as initial data. The mechanism of unidirectional block was not investigated. In this paper we simulate the FHN model on a ring with heterogeneities in the parameters and show that a rotating wave steady-state evolves from spatially homogeneous (trivial) initial data via a unidirectional block. In [4,8] propagation failure occurs in a spatial region described as a ‘‘gap” in which diffusion and/or excitability are altered. In addition, Ref. [13] reports that for a model of cardiac excitation on a ring with a spatial gap (at which there are heterogeneities in cellular coupling or size properties) a locally applied stimulus at the gap may result in a rotating wave on a ring. The main region of parameter heterogeneity in this paper will also be referred to as the gap. Here we are concerned with the transient behavior of the equations and, in particular, the parameters which result in unidirectional block. Our

E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.02.029


J.G. Alford / Applied Mathematics and Computation 216 (2010) 1337–1348

simulations show that a small asymmetry in the gradients of diffusion and stimulus at the border of the gap yields unidirectional propagation failure and a rotating wave steady-state. In order to understand one-way propagation failure at the border of the gap, the dynamics of the solutions to the equations near the gap must be understood. We do not analyze the FHN reaction–diffusion model here, but instead consider a system of four coupled cells with parameters that control unidirectional and bidirectional propagation. In [1,4,17] two and three coupled excitable cells are used to understand the dynamics of larger discrete and spatially continuous models. Unidirectional block and arrhythmia in cardiac cells has been studied both experimentally and with mathematical models using pairs of coupled cells [18–20]. Similarly, we consider a spatially discrete (four cell) model with a parameter to control gap dynamics that is motivated by our simulations on the ring. We use the McKean [21] (piece-wise constant) dynamics. The middle two cells make up the gap. A stimulus of sufficient magnitude is applied to the gap cells and they reach threshold from zero initial data. The outer two cells are coupled to the gap cells and are excitable (but do not recover). Unidirectional block occurs when one of the outer cells achieves threshold and the other does not. By considering a small time scale of recovery we neglect asymptotically small terms in the solutions. This yields a set of conditions and a relation between stimulus and coupling which result in unidirectional propagation. Arrhythmic behavior in cardiac cells is often (and naturally) considered to be a consequence of auto-oscillatory cell behavior (e.g. [11,20,22]). Throughout this paper we will only consider non-oscillatory kinetics as a mechanism for unidirectional block just as in the experimental work of [18] on pairs of rabbit ventricular cells. We show that the critical parameter region of stimulus over which unidirectional block may occur may be extremely narrow. Small perturbations of the stimulus away from this region results in no propagation (below criticality) or bidirectional propagation (above criticality). In Section 2 simulations of unidirectional block using the FitzHugh–Nagumo equations on a ring are shown. In Section 3 the four cell system with McKean kinetics is analyzed. 2. Simulations of unidirectional block on a ring The equations that will be considered in this section are a system of reaction–diffusion equations with FitzHugh–Nagumo dynamics. They are

  @v @ @v þ f ðv Þ  w þ IðhÞ; DðhÞ ¼ @t @h @h

@w ¼ b1 v  b2 w; @t


where f ðv Þ ¼ v ð1  v Þða  v Þ; a 2 ð0; 1=2Þ; p < h 6 p and t P 0. The geometry is one-dimensional circular so that (2.1) will be subject to 2p-periodic boundary conditions of the form

v ðpÞ ¼ v ðpÞ; v h ðpÞ ¼ v h ðpÞ;

wðpÞ ¼ wðpÞ:

The initial data will be v ðh; tÞ  0 and wðh; tÞ  0. The parameter functions for diffusion DðhÞ and stimulus IðhÞ are continuous, piecewise linear of the following form. Let 0 < a1 < b1 < p; 0 < a2 < b2 < p and

IðhÞ ¼ I1 ; IðhÞ ¼ I2 ;

h 2 ðp; b2 Þ [ ðb1 ; pÞ; h 2 ða2 ; a1 Þ;

DðhÞ ¼ d1 ; DðhÞ ¼ d2 ;




0.7 500










0.2 200

D 0.01

0.1 0

100 0



-0.1 -0.2 −3







0 −4









Fig. 1. The left-hand plot is a surface plot depicting v ðh; tÞ from solutions of (2.1) where the horizontal axis is h and the vertical axis is t. The bar to the right of this plot shows the magnitude of v. The corresponding parameter functions are shown in the right-hand plot. The horizontal axis is h and the vertical axis is either I or D as indicated in the plot. The parameter functions have values as in (2.3) and (2.4).


J.G. Alford / Applied Mathematics and Computation 216 (2010) 1337–1348


0.8 0.7


0.6 500

0.5 0.4






0.3 300

0.2 0.1


0 100 0




-0.1 -0.2 −3






0 −4










Fig. 2. Just as in Fig. 1 except the parameter functions (shown on the right) have the values from (2.5) and (2.6).


0.7 0.6


600 0.8 500

0.5 400

0.4 0.3


0.2 0.1








0 100 0

-0.1 -0.2 −3








100 0

-0.2 −3







Fig. 3. A surface plot depicting v ðh; tÞ from solutions of (2.1) with (2.2). The data is as in (2.3), (2.4) with I2 ¼ 0:02031 (left) and I2 ¼ 0:02033 (right). On the left, the stimulus is not large enough to initiate a wave. On the right, two rotating waves propagate in opposite directions and collide. In both cases, the steady-state is a spatially inhomogeneous stationary solution.

where I1 < I2 and d1 – d2 are small and positive. At the borders of the gap in the intervals ðb2 ; a2 Þ and ða1 ; b1 Þ, lines connect I1 ; I2 and d1 ; d2 . Figs. 1 and 2 show examples of these parameter functions (on the right) and the results of a simulation of (2.1) where for Fig. 1

a2 ¼ 0:34056; I1 ¼ 0;

b2 ¼ 0:55;

d1 ¼ 0:008;

a1 ¼ 0:39528;

I2 ¼ 0:020320310;

b1 ¼ 0:5;


d2 ¼ 0:005;


and for Fig. 2

a2 ¼ 0:20; I1 ¼ 0;

b2 ¼ 0:55;

d1 ¼ 0:0025;

a1 ¼ 0:3250;

b1 ¼ 0:5;

I2 ¼ 0:02157890625;

d2 ¼ 0:008:

ð2:5Þ ð2:6Þ

The method of lines was used to compute the solutions of the partial differential Eq. (2.1) shown in Figs. 1 and 3. The space variable h was discretized with an equally spaced grid of size Dh ¼ 2p=241  0:026 and MATLAB’s ode23s integration function was used to integrate in time along each hn . Standard finite difference schemes were used to approximate v h and v hh in (2.1). Numerical investigations of (2.1) with (2.2) showed that if all parameters but I2 were fixed, unidirectional rotating waves resulted only if I2 was in a very small interval. For example, with data as in (2.3) and (2.4) unidirectional rotating waves resulted only if I2 obeyed 0:02031 < I2 < 0:02033. If I2 < 0:02031 the steady-state was stationary, sub-threshold, non-constant as in Fig. 3(left). Conversely, if I2 > 0:02033 the steady-state was a bi-directionally propagating wave shown in Fig. 3(right). The critical values of I2 were determined by binary search.


J.G. Alford / Applied Mathematics and Computation 216 (2010) 1337–1348

3. Coupled cell model We will now analyze a discrete cell model of unidirectional block. For the continuous (diffusive) system, unidirectional block was caused by spatial asymmetries in diffusion and stimulus as depicted in Figs. 1 and 2. The mechanisms which lead to the block are dynamical and the steady-state is a rotating wave. In this section we consider the following semi-coupled system of linear ordinary differential equations representing four excitable ‘‘cells”

u01 ¼ dðu2  u1 Þ þ gðu1 Þ;


u02 ¼ d12 ðu1  u2 Þ þ d32 ðu3  u2 Þ þ gðu2 Þ  w2 þ I;

w02 ¼ bu2 ;

u03 ¼ d23 ðu2  u3 Þ þ d43 ðu4  u3 Þ þ c½gðu3 Þ  w3 þ I; u04 ¼ lðu3  u4 Þ þ gðu4 Þ;


w03 ¼ bu3 ;

ð3:3Þ ð3:4Þ

where b; d12 ; d32 ; d23 ; d43 ; l, and c are positive constants and the dynamics are as in the McKeon model where gðui Þ ¼ Hðui  aÞ  ui is the Heaviside function

Hðui  aÞ ¼

0; ui < a; 1; ui P a:

Note that the middle two Eqs. (3.2) and (3.3) include recovery variables w2 and w3 while the first and last equations do not recover. For convenience, each of (3.1)–(3.4) will be referred to as a cell. This discrete model is motivated in part by the simulations of the diffusive systems examples of which are shown in Figs. 1 and 2. If we consider the diffusion model to consist of a large number of coupled cells, then cells in the diffusive gap interval b2 < x < b1 are subjected to a stimulus and either weakly or strongly coupled to neighboring cells (depending on the value of the diffusion constant in the interval). The diffusive model has two intervals with non-zero gradients of diffusion and stimulus on the borders of the gap interval. Similarly, in the discrete model cell 3, whose dynamics are defined by (3.3), is a border zone cell which yields unidirectional block. Gradients in coupling strength in the diffusive gap interval are modeled in the discrete system by varying the coupling strengths in the forward coupling from cell 2 to cell 4 and cell 2 to cell 1. In order to analyze the discrete model we will consider a simplified system. If the coupling ratios obey

d12 =d32  1;

d43 =d23  1;


then the middle two cells are strongly coupled relative to the end cells just as in the parameter set (2.5) and (2.6) for the simulation shown in Fig. 2. The time courses shown on the left in Fig. 4 shows an example of unidirectional block for the discrete model (3.1)–(3.4) when its parameter set obeys (3.5). In this case, cell 1 reaches threshold and cell 4 does not. After numerically solving the system for a range of parameter values I and d (all others being fixed with l ¼ 0:9 d), the critical parameter values for which unidirectional block occurs is depicted (circles) on the right in Fig. 4. Here b ¼ 0:1; c ¼ 0:2; d32 ¼ d23 ¼ 0:003; d12 ¼ d43 ¼ 0:0003. For the same values of b; c, and d23 , Fig. 4 on the right also shows critical parameter value estimates (stars) for

d12 ¼ d32 ¼ d43 ¼ 0:





6.5 6




5 2




0 −1 0

3.5 5










3 0.02







Fig. 4. On the left are time courses which represent unidirectional block for the system (3.1)–(3.4) and a parameter set that obeys (3.5). On the right is a parameter space plot showing a region (UB) in parameter space I vs. d (all other parameters fixed with l ¼ 0:9d) of unidirectional block. The circles are from (3.5) and the stars from (3.6).

J.G. Alford / Applied Mathematics and Computation 216 (2010) 1337–1348


The critical parameter region for I and d from (3.6) approximates that from (3.5) when the region of block (UB) for a fixed value of d is large relative to the difference in the predictions of I from (3.5) and (3.6) as in Fig. 4. We will analyze the following partially coupled system which is essentially system (3.1)–(3.4) with (3.6)

u01 ¼ dðu2  u1 Þ þ gðu1 Þ;


u02 ¼ gðu2 Þ  w2 þ I;


w02 ¼ bu2 ;

u03 ¼ mðu2  u3 Þ þ c½gðu3 Þ  w3 þ I; u04

w03 ¼ bu3 ;

¼ lðu3  u4 Þ þ gðu4 Þ;

ð3:9Þ ð3:10Þ

where the coefficient m ¼ dð1  cÞ. The coupling in cells 2 and 3 is motivated by the diffusive system’s parameter functions described in (2.2). If c ¼ 1 and l ¼ d the system is symmetric: time courses for the pairs of cells 1, 2 and 3, 4 will be identical and there will be bidirectional propagation just as the diffusive equations yield bi-directionally propagating waves when b1 ¼ b2 and a1 ¼ a2 . If c – 1, we expect that certain parameter values may yield a unidirectionally propagating signal as the parameter set is asymmetric in this case. In the remainder of this section we will analyze the dynamics of (3.7)–(3.10) in order to determine parameter values d; l; I, and c for which unidirectional propagation may occur. The following assumptions on the parameters will be used 2

b < d ð1 þ dÞ1 ;

0 < a < 1=4;

b=a  d < a=2;


a < c < ð1  d  2bÞð1  dÞ ; 1 1 2 < I < ðd  3Þ; 2

ð3:11Þ ð3:12Þ ð3:13Þ

Cells 2 and 3 have unique equilibria ðu2 ; w2 Þ ¼ ðu3 ; w3 Þ ¼ ð0; IÞ while the equilibria for cells 1 and 4 can be the sub- or suprathreshold values

ui ¼ 0 () ui < a and ui ¼ ð1 þ dÞ1 () ui P a; for i ¼ 1; 4. Unidirectional propagation occurs when u1 ð1Þ ¼ ð1 þ dÞ1 and u4 ð1Þ ¼ 0. The dynamics of cells 2 and 3 will be considered in Sections 3.1 and 3.3. These results will be used to estimate the maxima for cells 1 and 4 in Sections 3.2 and 3.4 respectively. The dynamical behavior is then used to find parameters which yield unidirectional block in Section 3.5. 3.1. Cell 2 dynamics For arbitrary initial data u2 ðt0 Þ ¼ u02 and w2 ðt0 Þ ¼ w02 the time course for u2 ðtÞ is the solution of (3.8) given by

u2 ðtÞ ¼ a1 ejn1 jðtt0 Þ þ a2 ejn2 jðtt0 Þ ; 1 jn1 jðtt 0 Þ 1 n1 e


1 jn2 jðtt 0 Þ 2 n2 e

w2 ðtÞ ¼ ba þ ba þ I þ H2 ðu2  aÞ; pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x ¼ 1  4b; n1 ¼ ð1 þ xÞ=2 < 0; n2 ¼ ð1  xÞ=2 < 0; 1

a1 ¼ x


0 bu2 n1 2

a2 ¼ x1




þ I þ H2 ðu2  aÞ ;

ð3:15Þ ð3:16Þ ð3:17Þ

i 0 0 bu2 n1 1  w2 þ I þ H 2 ðu2  aÞ :


Since b  1, both n1 ; n2 are real-valued and u2 is non-oscillatory. If u2 is initially sub-threshold, a large enough stimulus parameter I can initiate a threshold response in the system. If u2 ð0Þ ¼ w2 ð0Þ ¼ 0 and uðtÞ < a, (3.14) becomes

u2 ðtÞ ¼ x1 Iðejn1 jt  ejn2 jt Þ;

t P 0:

ð3:19Þ 2

A Taylor expansion of x is that x ¼ 1  2b  2b     ; from which 2


jn1 j ¼ b þ b þ    ;

jn2 j ¼ 1  b  b     ;


1  4b < x < 1  2b;

b < jn1 j < 2b;



1  2b < jn2 j < 1  b:

These bounds yield a lower function, denoted f2 ðtÞ, where u2 ðtÞ > f2 ðtÞ and

f2 ðtÞ ¼ ð1  2bÞ1 Iðe2bt  eð12bÞt Þ:   1 This function is maximum at t M lnðg2 Þ and f2 t M ¼ Ir where 2 ¼ ð1  4bÞ 2

r ¼ ð1  2bÞ1 g22b=ð14bÞ  g2ð12bÞ=ð14bÞ ; g2 ¼ 2b=ð1  2bÞ: A numerical investigation shows that for 0 < b  1 the constant r decreases with b and r ¼ 1 in the limit as b # 0. A suf  ficient condition that u2 ðtÞ reaches threshold is that f2 t M > a which follows if I > ar1 > a which will be assumed. 2


J.G. Alford / Applied Mathematics and Computation 216 (2010) 1337–1348

In order to estimate a time to threshold for cell 2, notice that t M 2  lnð2bÞ and the expansion of n1 in (3.20) shows that t  jn1 j1 . Thus for t 6 OðtM tM 2 2 Þ; u2 from (3.19) is approximated by v2 ðtÞ ¼ Ið1  e Þ. After solving v2 ðtÞ ¼ a, an approximation for the time to threshold for u2 ðtÞ, denoted ^t 2 , is

^t2   lnð1  a=IÞ:


The assumptions on a and I in (3.11) and (3.13) shows that the time to threshold for cell 2 with trivial initial data is small, ^ as 1  a=I > 7=8 and ^t2 <  lnð7=8Þ  0:13. In this case, et2  1  ^t 2 ; v2 ð^t 2 Þ  I^t 2 , and since w02  bv02 ðtÞ, which yields that w02 ðtÞ  bv02 ðtÞ and w2 ð^t 2 Þ  bI^t 22 =2  a using assumptions on d and I in (3.11) and (3.13). To simplify the analysis when I is large and ^t2 small, the initial conditions will be taken as

u2 ð0Þ ¼ a;

u1 ð0Þ ¼ u3 ð0Þ ¼ u4 ð0Þ ¼ w2 ð0Þ ¼ w3 ð0Þ ¼ 0:

^ 2 ðtÞ Substitute u02 ¼ a; H2 ðu2  aÞ ¼ 1, and w02 ¼ 0 in (3.14), (3.17), and (3.18) and denote u2 ðtÞ as u

^ 1 ejn1 jt þ a ^ 2 ejn2 jt ; ^ 2 ðtÞ ¼ a t P 0; u

1 1 a^ 1 ¼ x ban2 þ I þ 1 ; a^ 2 ¼ x1 ban1 1 þIþ1 :

ð3:23Þ ð3:24Þ

3.2. Cell 1 dynamics ^ 2 ðtÞ in (3.23) and The sub-threshold dynamics for cell 1 will be determined by substituting the expression for u Hðu1  aÞ ¼ 0 into (3.7). For initial data u1 ð0Þ ¼ 0, the time course for cell 1 is then

jn1 jt jn2 jt ^ 1 c1 ^ 2 c1 ; u1 ðtÞ ¼ ^c1 eð1þdÞt þ d a þa 11 e 12 e

t P 0;



  ^ 1 c1 ^ 1 ^c1 ¼ d a 11 þ a2 c12 ;

c11 ¼ 1 þ d  jn1 j; c12 ¼ 1 þ d  jn2 j:


Now the expansions (3.20) with (3.24) yield that

a^ 1 ¼ ð1 þ IÞð1 þ 2bÞ þ Oða bÞ; a^ 2 ¼ ð1 þ I  aÞð1 þ 2bÞ þ Oða bÞ:


Note that (3.20) shows c11  1 þ d and c12  d þ b so that

^c1  ð1 þ I  aÞð1 þ b=dÞ1  dð1 þ dÞ1 ð1 þ IÞ:


The assumptions (3.11) and (3.13) show that ^c1 is positive and u1 ðtÞ > f1 ðtÞ where

jn1 jt jn2 jt ^ 1 c1 ^ 2 c1 f1 ðtÞ ¼ d a þa ; 11 e 12 e

t P 0:


1 The maximum of f1 ðtÞ is at tM ln W1 1 and 1 ¼ x

h i   jn1 j=x jn2 j=x ^ 2 c1 ^ 1 c1 f1 t M þa ; ¼d a 1 11 W1 12 W1

W1 ¼ 

^ 1 c12 jn1 ja : ^ 2 c11 jn2 ja 


A sufficient condition so that u1 ðtÞ reaches threshold is that f1 t1


bd ð1 þ IÞ 12 < bd; ð1 þ I  aÞð1 þ dÞ 11


  12 >  ln tM bd ; 1 11


P a. Use the expansion (3.20) with (3.27) to get



after neglecting terms of order b ð1 þ IÞ and using the bound a=ð1 þ IÞ < 1=12 which follows from assumptions (3.11) and (3.13). This then shows that M

ð1 þ IÞet1 < 12 bd ð1 þ IÞ=11  a;


    M and (3.28) and (3.32) shows ^c1 eð1þdÞt1  a. Thus (3.29) approximates u1 ðtÞ < a for t ¼ O tM and f1 t M approximates max1 1 imum sub-threshold u1 ðtÞ.  M In order to approximate f1 t 1 , we will derive bounds on the various constants in (3.30). First use the approximation in (3.31) and assumptions (3.11) and (3.13) to see that

bd ð1 þ dÞ1 < W1 < bd ð1  a=3Þ1 ð1 þ dÞ1 < 1:


2 j=x It follows that Wjn < W1 as jn2 j > x from (3.21). Furthermore (3.21) shows that jn1 j  b; x  1 and from (3.27) we get 1       ^ ^ a1  1 þ I; a2 < 0. Thus (3.30) and the bounds in (3.33) yields f1 tM1 < f1 tM1 < f1 tM1 where

  ¼ dð1 þ dÞ1 ð1 þ IÞv1 ðdÞ; f1 t M 1   ¼ dð1 þ dÞ1 ð1 þ IÞv1 ðdÞ; f1 t M 1

ð3:34Þ ð3:35Þ

J.G. Alford / Applied Mathematics and Computation 216 (2010) 1337–1348



v1 ðdÞ ¼

b bd ; 1þd

v1 ðdÞ ¼

bd ð1  a=3Þð1 þ dÞ

b ð3:36Þ


after ignoring terms of order bd ð1 þ IÞ or less. 3.3. Cell 3 dynamics The sub-threshold solution of the non-homogeneous initial value problem for cell 3 using (3.9) with u3 ð0Þ ¼ 0 and ^ 2 from (3.23) and Hðu3  aÞ ¼ 0. It is given by w3 ð0Þ ¼ 0 is obtained after substituting u

u3 ðtÞ ¼ h1 ejk1 jt þ h2 ejk2 jt þ c11 ejn1 jt þ c12 ejn2 jt ;


jk1 jt jk2 jt w3 ðtÞ ¼ bk1 þ bk1 þ c21 ejn1 jt þ c22 ejn2 jt þ I; 1 h1 e 2 h2 e



qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi /2  4cb;

/ ¼ c þ m;

k1 ¼ ð/ þ qÞ=2 < 0;


k2 ¼ ð/  qÞ=2 < 0;


^ 1 ð/  jn1 j  cb=jn1 jÞ1 ; c21 ¼ bjn1 j1 c11 ; c11 ¼ ma


^ 2 ð/  jn2 j  cb=jn2 jÞ1 ; c22 ¼ bjn2 j1 c12 ; c12 ¼ ma and

h i h1 ¼ cq1 I þ b jk2 j1 ðc11 þ c12 Þ þ c21 þ c22 ; h i h2 ¼ cq1 I þ b jk1 j1 ðc11 þ c12 Þ þ c21 þ c22 :


Note that (3.11) and (3.12) imply /2 =c > 4b if a > 4b as / > c > a. Thus it will be assumed that a > 4b, q is real-valued, and the time courses are non-oscillatory. It is seen from (3.16) and (3.40) that jk1 j < jn1 j precisely when x þ /  q < 1. Differentiation shows

d ðx þ /  qÞ ¼ 2cq1  2x1 < 0 db After substituting q ¼


q2 c2 > x2 :


qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi /2  4cb and x ¼ 1  4b and noting that /2 > c2 þ 2mc, then (3.43) holds if d > 2b which is as-

sumed. Therefore x þ /  q decreases with b and must be less than 1 (as x þ /  q ¼ 1 when b ¼ 0) and jk1 j < jn1 j. Next, (3.16) and (3.40) show that jk2 j < jn2 j if and only if / þ q < 1 þ x which follows from the bound on x in (3.21) and assumed bound on c in (3.12). Finally, use (3.11) and (3.21) to see that jn1 j < 2b < a < / shows jn1 j  jk2 j. In summary, the eigenvalues which define the time constants in (3.37) compare as

jk1 j < jn1 j  jk2 j < jn2 j:


In order to approximate the various constants in (3.37), first Taylor expand q to get

q ¼ /ð1  2cb/2  2c2 b2 /4    Þ ¼ /  2cb/1  Oð2b2 /1 Þ; 2 2 jk1 j ¼ cb/1 þ Oðb /1 Þ; jk2 j ¼ /  cb/1  Oðb /1 Þ:

ð3:45Þ ð3:46Þ


From (3.20) it is seen that after neglecting terms Oðb Þ; cb=jn1 j  cð1  bÞ and cb=jn2 j  cbð1  bÞ1 , so that

/  jn1 j  cb=jn1 j  mð1  b=dÞ;

/  jn2 j  cb=jn2 j  /  1;

and (3.27), (3.41), (3.47) with the identity m=ð1  /Þ ¼ dð1  dÞ

c11 ¼ ð1 þ IÞð1 þ 2bÞð1  b=dÞ1 þ OðbÞ; c12 ¼ dð1 þ I  aÞð1 þ 2bÞð1  dÞ





c21 ¼ c11 ð1  OðbÞÞ;

þ OðbÞ; c22 ¼ bc12 ð1 þ OðbÞÞ:


The assumptions (3.11) and (3.12) show b  a < c < / and

bð1 þ IÞð1  b=dÞ1 

a a dð1 þ IÞð1  b=dÞ1 < ð1  b=dÞ1 ; 2 2

and dð1  dÞ1 < 1=7 so that (3.42), (3.46) and (3.48) give the approximations



J.G. Alford / Applied Mathematics and Computation 216 (2010) 1337–1348

q  /; bjk1 j1  /c1 ; bjk2 j1  b/1 ; 1

ð3:50Þ 1

c11  ð1 þ IÞð1  b=dÞ ; c12  dð1 þ I  aÞð1  dÞ ;

h1  c/1 I  c11 ð1  b/1 Þ ; h2  c/1 I þ c11 ½c/1  1  c12 :

ð3:51Þ ð3:52Þ

The approximations in (3.52) give that

h1 þ c11  c11 ð1  c/1 þ cb/2 Þ þ c/1 I; h2  h1  c11  c12  cb/2 ;


and the assumptions (3.11) and (3.13) with (3.51) and (3.53) yield

1 þ I  c/1 < h1 þ c11 ;

0 < c12 < 1=2:


The approximations above can be used to derive bounds for sub-threshold u3 ðtÞ in (3.37) as follows. Firstly, h1 þ c11 > 0 and h2 < 0 (as c < /) with (3.44) shows

h1 ejk1 jt þ c11 ejn1 jt < ðh1 þ c11 Þejk1 jt < h1 þ c11 :



u3 ðtÞ ¼ h1 þ c11 ejn1 jt þ h2 ejk2 jt ;

u3 ðtÞ ¼ h1 þ c11 þ h2 ejk2 jt þ c12 ejn2 jt ;


and inspection of (3.37), (3.44) and (3.55) shows u3 ðtÞ 6 u3 ðtÞ 6 u3 ðtÞ After Taylor expanding the exponentials ejn1 jt and ejk2 jt , and using h2 < 0 it is seen that u3 ðtÞ > f3 ðtÞ where

f3 ðtÞ ¼ h1 þ c11 ð1  jn1 jtÞ þ h2 ð1  jk2 jt þ jk2 j2 t2 =2Þ: Now (3.44), (3.51) and (3.52) show jh2 j > c/1 I and c11 jn1 j  jh2 j jk2 j from which

 maxff3 g  f3 jk2 j1  h1 þ c11 þ h2 =2  ðh1 þ c11  c12 Þ=2:  Then maxff3 g  f3 jk2 j1 > a, which follows from assumptions (3.11) and (3.13), and the bound (3.54). Thus u3 ðtÞ P f3 ðtÞ and u3 ðtÞ ¼ a for some t < jk2 j1 . Denote ^t 3 to be the time for which cell 3 reaches threshold so that u3 ð^t 3 Þ ¼ a and ^t 3 < jk2 j1  jn1 j1 . For t ¼ Oð^t3 Þ, the lower bound for sub-threshold u3 ðtÞ in (3.56) yields u3 ðtÞ ¼ a if

a  h1 þ c11 þ h2 ejk2 jt


^t 3 < jk2 j1 ln h1 þ c11 þ c12 ; h1 þ c11  a


using the approximation in (3.53) for h2 . Now use (3.46) to approximate jk2 j  / > c and the bounds (3.54) and (3.57) show that

^t3 < /1 ln h1 þ c11 þ c12 < c1 ln I þ 1=2 < 1; Ia h1 þ c11  a


where the final inequality follows from the assumptions (3.11) and (3.13) which show that I > 2 > ð1=2 þ aec Þðec  1Þ1 . The estimates (3.44) and (3.58) and comparing sub-threshold u3 ðtÞ from (3.37) with the upper bound for sub-threshold u3 ðtÞ in (3.56) shows that u3 ðtÞ  u3 ðtÞ for t ¼ Oð^t3 Þ. Furthermore, jn2 j < 1 and jk2 j < 1 show that ejn2 jt  1  jn2 jt and ejk2 jt  1  jk2 jt. Substitute these linearizations into u3 ðtÞ and equate the result to a to get that

^t3  ðh1 þ c11 þ h2 þ c12  aÞðc12 jn2 j þ h2 jk2 jÞ1 : Next substitute jk2 j  /; jn2 j  1 and use (3.51)–(3.53) and the identities 1  / ¼ ð1  dÞð1  cÞ and see that

c12 jn2 j þ h2 jk2 j  c12 þ /h2  c12 ð1  /Þ  c11 ð/  cÞ  cI < cI;


m ¼ /  c ¼ dð1  cÞ to ð3:60Þ

as c11 > 1 þ I and c12 < dð1  dÞ1 ð1 þ IÞ. Finally, combine (3.53), (3.59) and (3.60) and assumptions (3.11) and (3.13) to see that

^t3  h1 þ c11 þ h2 þ c12  a < a < 1 : c12 jn2 j þ h2 jk2 j cI 2


The approximations (3.53) and (3.54) with sub-threshold u3 ðtÞ  u3 ðtÞ from (3.56) and assumptions (3.11) and (3.13) M ^ 2 jet1  a. From these observations, the approxshows that jw3 ð^t3 Þj  jbu03 ð^t3 Þj  a. Inspection of (3.27) and (3.32) shows ja   jn jt ^ 1 e 1 will be used to compute u3 ðtÞ for t ¼ O tM for the non-homogeneous initial value imations w3 ð^t3 Þ  0 and u2 ðtÞ  a 1 problem (3.9) with Hðu3  aÞ ¼ 1; u3 ð^t3 Þ ¼ a. This gives an approximation of suprathreshold u3 ðtÞ when c – 0 as

J.G. Alford / Applied Mathematics and Computation 216 (2010) 1337–1348

  ^ 3 ðtÞ ¼ ^h1 ejk1 jðt^t3 Þ þ ^h2 ejk2 jðt^t3 Þ þ c11 ejn1 jt ; t ¼ O tM u3 ðtÞ  u 1 ; h  i ^h1 ¼ cq1 1 þ I  bjn1 j1 c11 ejn1 j^t3  bjk2 j1 a  c11 ejn1 j^t3 ; h  i ^h2 ¼ cq1 1 þ I  bjn1 j1 c11 ejn1 j^t3  bjk1 j1 a  c11 ejn1 j^t3 :


ð3:62Þ ð3:63Þ

The analysis of cell 3 shows how c controls certain asymmetries in the time courses. The estimate (3.22) of time to threshold for cell 2 is that

^t2 ¼  lnð1  a=IÞ ¼ a=I þ Oða2 =I2 Þ; and (3.61) then shows that ^t 2  c^t 3 . There is an effective delay in the time to threshold for cell 3 (relative to cell 2) induced by c when 0 < c < 1. In addition, the time length over which cell 3 remains suprathreshold is controlled by c. This can be seen in the single cell dynamics (set d ¼ 0 in (3.9)) after re-scaling the time variable as s ¼ ct to get

du3 =ds ¼ u3 þ Hðu3  aÞ  w3 þ I;

dw3 =ds ¼ c1 bu3 :

The suprathreshold decay rate (fast time scale) is increased when 0 < c < 1 as can be seen from the suprathreshold single cell dynamics computed in (3.23) with time constants as in (3.20). In summary, when 0 < c < 1 there will be shorter time intervals of excitability and longer time intervals of recovery for cell 3 relative to cell 2 which is apparent in the time courses in Fig. 4. 3.4. Cell 4 dynamics The sub-threshold dynamics for u4 ðtÞ may be determined with u3 ðtÞ as in (3.37) and (3.62). Firstly, if the approximation u3 ðtÞ  u3 ðtÞ from (3.56) is substituted into (3.10) with u4 ð0Þ ¼ 0 and Hðu4  aÞ ¼ 0, the result is that

h1 þ c11 h2 ejk2 jt c12 ejn2 jt ; u4 ðtÞ  g4 eð1þlÞt þ l þ þ 1þl 1 þ l  jk2 j 1 þ l  jn2 j

0 6 t 6 ^t 3 ;

where g4 is defined so that the right side is zero when t ¼ 0. To approximate u4 ð^t3 Þ, use the estimate ^t3 from (3.61) and sub^ stitute ert3  1  r^t3 þ ðr^t 3 Þ2 =2 with r ¼ jn2 j  1; r ¼ jk2 j  /, and r ¼ 1 þ l gives

u4 ð^t 3 Þ  l^t3 ½h1 þ c11 þ h2 þ c12  Oð^t3 Þ  0;


using the estimates (3.53) and (3.61). Next, let Hðu4  aÞ ¼ 0 and u3 ðtÞ be from (3.62) when t > ^t 3 and u4 ðtÞ < a. After substituting into (3.10) and solving with (3.64), the result is that ð1þlÞðt^t 3 Þ

u4 ðtÞ  ^c4 e


" c11


jn1 jt





jk1 jðt^t 3 Þ





# jk2 jðt^t 3 Þ





h i jn1 j^t3 ^ 1 ; ^c4 ¼ l c11 u1 þ ^h1 u1 3 e 1 þ h2 u2

u1 ¼ 1 þ l  jk1 j; u2 ¼ 1 þ l  jk2 j; u3 ¼ 1 þ l  jn1 j;


  valid for t ¼ O tM 1 . Estimates of the various constants in (3.65) will now be determined. Substitute the approximations from (3.50) into (3.63) to see that

^h1  cq1 1 þ I  c11 ð1 þ b/1 Þ ;

^h2  cq1 1 þ I þ c11 ð/c1  1Þ  a/c1 ;

ð3:67Þ ð3:68Þ

after neglecting the term a b/1 < b (as / > a from (3.11) and (3.13)). Now use jn1 j  b; jk1 j < jn1 j, and jk2 j  / from (3.20), (3.44) and (3.46) to see that

u1  1 þ l; u2  1 þ l  /; u3  1 þ l:


Now it is easily seen from assumptions (3.11), (3.13) and (3.48) that ^ h1 < 0 and ^ h2 < 0 as / > c. Thus after substituting (3.67)–(3.69), we get that 1 1 ^ 1 ^ 1 ^ ^ c11 u1 < 1=2; 3 þ h1 u1 þ h2 u2 < ð1 þ lÞ ðc 11 þ h1 þ h2 Þ < a/c 1 1 1 1 1 1 ^ ^ c11 u3 þ h1 u1 þ h2 u2 > c11 q ð1 þ lÞ ðm  cb/ Þ  ð1 þ IÞ > ð1 þ IÞ;

where we have used m ¼ /  c and the assumptions (3.11) and (3.13) show that a < c < / < 1=2; a d > b, and c < q. These bounds with (3.32) show that


J.G. Alford / Applied Mathematics and Computation 216 (2010) 1337–1348 ^ j^c4 jeð1þlÞðtt3 Þ  a;


as ^t 3 < 1=2  tM 1 from the estimates (3.31) and (3.61). 1 ^ Moreover, (3.46) shows that jk1 j < b, tM 1  t 3  jk1 j , and (3.65) and (3.70) yield the estimate

h i 1 jn1 jt jk2 jt u4 ðtÞ  f4 ðtÞ ¼ l ^h1 u1 þ ^h2 u1 ; 1 þ c 11 u3 e 2 e


1 1 which is maximum at tM 4 ¼ b4 ln W4 ,

h i jk2 j=b4 1 jn1 j=b4 ^ 1 ; f4 ðtM þ ^h2 u1 4 Þ ¼ l h1 u1 þ c 11 u3 W4 2 W4 b4 ¼ jk2 j  jn1 j;


W4 ¼ ðjn1 jc11 u2 Þ=ðjk2 j^h2 u3 Þ:

We will next derive bounds on (3.72). First use jn1 j  b which follows from (3.20).Then b=/ < b=d  0, which follows from (3.11), and (3.46) shows that jk2 j  /. Then (3.68), (3.69) and (3.72) yield the approximation

bð1 þ l  /Þ : 1  1 ð1 þ lÞ c ð1 þ I  a/c1 Þc1 11 þ /c



Then from (3.51) and assumptions (3.11) and (3.13) we use c11  1 þ I < 3 to get bounds


a/ < ð1 þ I  a/c1 Þc1 11 < 1; 3c

so that (3.73) yields that

bð1 þ l  /Þ b ð 1 þ l  /Þ < W4 < : /ð1 þ lÞ /ð1  a=3Þð1 þ lÞ jk j=b4

Now note that W4 2

^h2 u1 W 2

jk2 j=b4 4




 W4 and from (3.72) it is seen that

1 1  ^h2 u1 2 W4  b/ ð1 þ lÞ ð1 þ IÞ  0;


using the previous approximations for jn1 j; c11 ; jk2 j, and u3 . Finally, using (3.72), (3.74), and (3.75) we get the bounds       f4 tM < f4 tM < f4 t M where 4 4 4

  ¼ lð1 þ lÞ1 ð1 þ IÞv4 ðd; c; lÞ; f4 t M 4  M f4 t 4 ¼ lð1 þ lÞ1 ð1 þ IÞv4 ðd; c; lÞ;

ð3:76Þ ð3:77Þ




















Fig. 5. The critical curve I ¼ j1 ðdÞ from (3.79) is the lower bound solid curve. The critical curve I ¼ j4 ðd; c; lÞ from (3.79) is the middle bound dashed curve when c ¼ 0:9 and the upper bound solid curve when c ¼ 0:16. Unidirectional block occurs when the parameters lie between the lower solid curve and the upper dashed curve when c ¼ 0:9 and the lower solid curve and the upper solid curve when c ¼ 0:16.

J.G. Alford / Applied Mathematics and Computation 216 (2010) 1337–1348



v4 ðd; c; lÞ ¼

 b=/ bð1 þ l  /Þ ; /ð1 þ lÞ

v4 ðd; c; lÞ ¼

bð1 þ l  /Þ /ð1  a=3Þð1 þ lÞ

b=/ ;



after ignoring terms of order b l d ð1 þ IÞ or less. 3.5. Parameter space estimates In this section we approximate the regions in parameter space for which unidirectional block will occur. Setting the lower bound for cell 1 in (3.34) and the upper bound for cell 4 in (3.77) to threshold a and solving for I gives the functions

 1 a 1þd    1; j1 ðdÞ ¼ f1 t M 1

 1 a 1þd    1; j4 ðd; c; lÞ ¼ f4 tM 4


as an estimate of the maximum critical I value for cell 1 and minimum critical I value for cell 4 to achieve threshold for a given set of parameters d, c, and l. Then the region in parameter space defined by

j1 ðdÞ < I < j4 ðd; c; lÞ


approximates the region for which there will be unidirectional block. The estimate (3.80) holds under the assumptions (3.11)–(3.13). We will also require that the errors in the bound for max M 1 imum cell 1, Oðb d ð1 þ IÞÞ, and maximum cell 4, Oðb l d ð1 þ IÞÞ, are small relative to threshold value a and that tM 4 ¼ O t1 . 1 1 1 1 ln W1 ; x  1, and (3.33) with t M Use tM 1 ¼ x 4 ¼ b4 ln W4 ; b4  /, and (3.74) to see that


   ð1  a=3Þð1 þ dÞ 1þd < tM ; 1 < ln bd bd



/1 ln

    /ð1  a=3Þð1 þ lÞ /ð1 þ lÞ 1 < / ln < tM : 4 bð1 þ l  /Þ b ð 1 þ l  /Þ


 M These bounds will be used to check that t M 4 ¼ O t1 . Fig. 5 depicts the region (3.80) for 0:045 < d < 0:05 when a ¼ 0:15; b ¼ 0:0001, and c is either c ¼ 0:16 or c ¼ 0:9. Here  M M the coupling  M M  constant for cell 4 is l ¼ 0:99 d. TheM bounds  M (3.81) and (3.82) yield that max t4 =t1  3 when c ¼ 0:16 and max t4 =t1  1 when c ¼ 0:9. It follows that t 4 ¼ O t1 . Furthermore, the relative errors in the approximations for j1 1 and j4 given by b l d ð1 þ IÞ=a;  0:002 and b dð1 þ IÞ=a  0:002 are small. 3.6. Conclusion In Section 2 we investigated numerical simulations of unidirectional block in the FitzHugh–Nagumo reaction–diffusion equation on a ring. For spatially homogeneous (zero) initial data and spatially heterogeneous diffusion and stimulus the steady-state is a rotating wave. For the diffusion and stimulus functions described in (2.2), the rotating wave occurs for an extremely narrow range of peak stimulus I. The kinetics are non-oscillatory (that is, in the absence of diffusion, the steady-state is independent of time). In Section 3 a discrete model of four coupled excitable cells which exhibit unidirectional block was analyzed. The solutions were approximated assuming the parameter which controls recovery time was small and a coupling and stimulus range was determined for which unidirectional block may occur. An important distinction between the diffusive and discrete models here is the cell-to-cell coupling schema. The boundary conditions for the diffusive model were periodic and the geometry was a ring. The interaction of the oppositely traveling wavefronts (clockwise and counter-clockwise) will influence the resulting steady-state. An analogous discrete model would be periodically coupled (i.e. cells 1 and 4 include recovery and are mutually coupled) and the analogous steady-state is oscillatory dynamics and a self-propagating excitation through each cell. For the diffusive model the asymmetries in the excitation and recovery times of the two opposite traveling fronts initiated by the stimulus will influence the steady-state behavior. In particular, a rotating wave will be established only if the leading cells (those in the direction of travel of the front) are excitable as the front approaches. This dependence was observed in [13]. Thus the diameter of the ring and the speed of the wave should be considered in the analysis of the diffusive model. For the discrete model it was discussed at the end of Section 3.3 how c influences the excitation and recovery time in cell 3 which may inhibit or promote a self-propagating impulse in the periodically coupled discrete system. Thus c may be used in future work to investigate the influence of excitation and recovery times on the steady-states in the periodically coupled discrete model.


J.G. Alford / Applied Mathematics and Computation 216 (2010) 1337–1348

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

V. Booth, T. Erneux, Understanding propagation failure as a slow capture near a limit point, SIAM J. Appl. Math. 55 (1995) 1372–1389. H. Ikeda, R. Mimura, Wave-blocking phenomena in bistable reaction–diffusion systems, SIAM J. Appl. Math. 49 (1989) 515–538. J.P. Keener, Propagation and its failure in coupled systems of discrete excitable cells, SIAM J. Appl. Math. 47 (2000) 556–572. T.J. Lewis, J.P. Keener, Wave-block in excitable media due to regions of depressed excitability, SIAM J. Appl. Math. 61 (2000) 293–316. S. Lubkin, Uni-directional waves on rings: models for chiral preference of circumnutating plants, Bull. Math. Biol. 56 (1994) 795–810. J.P. Pauwelussen, Nerve impulse propagation in a branching nerve system: a simple model, Physica D 4 (1981) 67–88. A. Rabinovitch, I. Aviram, N. Gulko, E. Ovsyshcher, A model for propagation of action potentials in non-uniformly excitable media, J. Theor. Biol. 196 (1999) 141–154. J. Sneyd, J. Sherratt, On the propagation of calcium waves in inhomogeneous medium, SIAM J. Appl. Math. 57 (1997) 73–94. C. Antzelevitch, Basic mechanisms of reentrant arrhythmias, Curr. Opin. Cardiol. 16 (2001) 1–7. L. Glass, M.E. Josephson, Resetting and annihilation of reentrant abnormally rapid heartbeat, Phys. Rev. Lett. 75 (1995) 2059–2062. H. Gonzalez, Y. Nagai, G. Bub, L. Glass, A. Shrier, Reentrant waves in a ring of embryonic chick ventricular cells imaged with a Ca2+ dye, Biosystems 71 (2003) 71–80. H. Ito, L. Glass, Theory of reentrant excitation in a ring of cardiac tissue, Physica D 56 (1992) 84–106. W. Quan, Y. Rudy, Unidirectional block and reentry of cardiac excitation: a model study, Circ. Res. 66 (1990) 367–382. J.G. Alford, G. Auchmuty, Rotating wave solutions of the FitzHugh–Nagumo equations, J. Math. Biol. 53 (2006) 797–819. R. FitzHugh, Impulses and physiological states in theoretical models of nerve membrane, Biophys. J. 1 (1961) 445–466. J.S. Nagumo, S. Arimoto, S. Yoshizawa, An active pulse transmission line simulating nerve axon, Proc. IRE 50 (1962) 2061–2071. G.B. Ermentrout, J. Rinzel, Reflected waves in an inhomogeneous excitable medium, SIAM J. Appl. Math. 56 (1996) 1107–1128. R. Joyner, H. Sugiura, R. Tan, Unidirectional block between isolated rabbit ventricular cells coupled by a variable resistance, Biophys. J. 60 (1991) 1038– 1045. M. Landau, P. Lorente, Conduction block and chaotic dynamics in an asymmetrical model of coupled cardiac cells, J. Theor. Biol. 186 (1997) 93–105. B.E. Peercy, J.P. Keener, Coupled cell model of border zone arrhythmias, SIAM J. Appl. Dyn. Sys. 4 (2005) 679–710. H.P. McKean, Nagumo’s equation, Adv. Math. 4 (1970) 209–223. R. Aliev, A. Panfilov, Modeling of heart excitation patterns caused by a local inhomogeneity, J. Theor. Biol. 181 (1996) 33–40.