Temporal stability of solitary impulse solutions of a nerve equation

TEMPORAL STABILITY OF SOLITARY IMPULSE SOLUTIONS OF A NERVE EQUATION JOHN A. FEROE, Department of Mathematics, Vassar College, Poughkeepsie, New York 12601 U.S.A. ABSTRACT We study a differential equation that models nerve impulse transmission. The nonlinearity is simplified to be piecewise linear in order to allow explicit solution. In general, two solitary impulse solutions are exhibited. The temporal stability of these solutions is analyzed by a technique that identifies the number of unstable modes. These results extend the results of Rinzel and Keller (1973, Biophys. J. 13:1313) by showing that the slower unstable solution has only one unstable mode, and that the fast solution, as conjectured, has no unstable modes and is therefore stable.

INTRODUCTION

The best known and most important mathematical model for the dynamics of nerve impulse propagation is the partial differential equation of Hodgkin and Huxley (10). Essential qualitative features of their model are preserved in a simplified equation introduced by FitzHugh (8,9) and Nagumo et al. (13). The nonlinearities in both of these equations have prevented exact solution, a fact which led McKean (12) to propose a piecewise linear form of the FitzHugh-Nagumo equation that can be explicitly solved. As studied by Rinzel and Keller (14, 15), this piecewise linear model, which we continue to investigate here, has provided valuable insight into the more elaborate models. The question of the stability of propagating impulse solutions has always been central in the study of such equations. An impulse is assumed to be subjected to a perturbation, which at the initial time is bounded along the length of the axon. If the perturbation grows in time, the impulse is unstable, while if it decays, it is stable. Recent work by Evans, (3-6), summarized in ref. 7, has characterized stability in a general setting. The application of these new techniques to the piecewise linear model constitutes an important aspect of the study of the model. It illustrates as well the basic idea of the method, since the constructions in this setting have a particularly transparent form. We consider the equation Ov/dt = d2v/0x2 - v - w + H(v - a), awlat = by - dw,

(1)

for positive constants a, b, and d, where H is the Heaviside step function BIOPHYSICAL JOURNAL VOLUME 21 1978

103

H(v - a) =.

I if V > a

'.0if Pa

The function v(x, t) models the transmembrane voltage of a nerve axon at distance x and time t, while w(x, t) is an auxiliary variable often called the recovery variable (9). Here we set d to be a small positive value, but observe that the results of this paper are qualitatively identical for d = 0, the value used by Rinzel and Keller. A nonzero value of d serves to introduce resistance in the dynamics of w and causes the rest solution v w 0 to be exponentially stable (4). A traveling wave solution of Eq. 1 is given by (v(x, t), w(x, 1)) = (vj(z), wj(z)), where z = x + ct, c > 0. Such functions of z are solutions of the related ordinary differential equation

czC =

vc

-vc - wc + H(vc - a)

cw' = bv -dw.

(2)

Correspondingly, bounded solutions to Eq. 2 can be viewed as solutions to Eq. 1 which represent waves of fixed form traveling to the left with velocity c. Only for limited values of the parameters a, b, and d will Eq. 2 exhibit such bounded solutions. Rinzel and Keller found there to be at most two values of c (denoted cf and c,) for which Eq. 2 has a solution with vc(z) 0 and wj(z) - 0 as z -x o, and which has the characteristic rising-falling-recovery solitary impulse form. For other parameter values 4-

06 -I 04 0.2 ~~ ~~1020

-5

-0.2 0-I 0

I

\

01

0?2 a

FIGURE I

03

/

-0.4 FIGURE 2

FIGURE 1 The values of c for which Eq. 2 has a solitary impulse solution. Here b = 0.2, d = 0.05, and a > 0. The upper portion of the curve gives values of Cf, the lower portion gives values of c5, and c, occurs at the point with a vertical tangent line. FIGURE 2 The graphs of the solitary impulses v,(z) for velocity values c, (solid line) and Cf (dotted line). The parameter values used here are a = 0.24, b = 0.2, d = 0.05, Cf = 1.0779, and c, =0.4720. The values for z1 where vc(zl) = a are approximately 1.30 for the slow impulse and 6.27 for the fast impulse.

104

BIOPHYSICAL JOURNAL VOLUME 21 1978

there is either a unique solitary pulse solution (with velocity denoted as c,) or no such solution. Periodic solutions are possible for some values of c, but are not considered in this paper. In Fig. 1 we show for b = 0.2 and d = 0.05 the dependence of c on values of a. (See ref. 15 for a more complete figure.) In Fig. 2 we show for a = 0.24, b = 0.2 and d = 0.05 both pulse solutions, with the leading wavefront spatially located with vM(O) = a. Both figures are characteristic. This distribution of solutions seems to be common to all widely accepted nerve models (see 1, 2,9). Huxley (11) conjectured for the Hodgkin-Huxley model that since only the faster impulse is seen physiologically, it must be stable, while the smaller slower impulse must be unstable. Rinzel and Keller (15) were able to demonstrate the instability of the slower pulse solutions by locating the growth factor of an unstable mode. They further provide figures showing the dependence of that growth factor on the parameter values. Inasmuch as a demonstration of the conjectured stability of the faster pulse requires showing the nonexistence of an unstable growth mode, that was not possible using their techniques. We present here a demonstration of the stability of the faster impulse solution for certain parameter values, and a demonstration that the slow impulse solution has only one unstable mode. While the results are valid only for the parameter values for which computations were done, we are confident that the results are completely representative. STABILITY OF IMPULSE SOLUTIONS

In this section we outline the mathematical theory and present the results on the stability of impulse solutions to Eq. 1. The reader is referred to the appendix for technical details on the mathematics, and to refs. 7 and 15 for expanded discussions of the methods. A solitary pulse solution (va, w,) has only two z values, z0 and z1, for which v'(z) = a (see Fig. 2). As z does not appear explicitly in Eq. 2, any translate of a solution is also a solution and thus we are free to position (va, w,) anywhere on the z-axis by a choice of z0. We follow Rinzel and Keller and take z0 = 0. Since the Heaviside function only changes at 0 and zi, on any one of the three intervals, (-e , 0), (0, zi), (z1, +o), the impulse (va, wc) is a solution of a completely linear system, that is, the combination of (three) vectors with exponential coefficients. Only one of these vectors has an exponential coefficient with positive real part and so to satisfy the requirement that it be bounded, (va, wc) on the interval (-oo ,0) must be a multiple of that function, in particular the multiple for which vc(0) = a. In fact (va, wc) is uniquely determined for all z since the full solution requires that the separate solutions on each interval match at 0 and z1 . This restriction accounts for the fact that so few values of c give rise to impulse solutions as well as provides the method for locating such c values. In the moving coordinate system (z, t) a pulse solution (vj(z), wc(z)) is a time-independent solution of

dv/dt = 02v/0z2 - c dv/dz - v - w + H(v - a) dw/dt = - c w/Oz + bv - dw, J. A. FEROE Solitary Impulse Solutions of a Nerve Equation

(3) 105

and its stability depends on the properties of the linearization of Eq. 3 about (vC, w), which is av/ot = 82v/0z2 - c av/Ot - (1 - b(z)/vC(0) + b(z - z1)/v'(z1)) v - w (4) aw/at = - caw/Oz + bv - dw, where 6 is the Dirac delta function. In particular if Eq. 4 has no solution of the form

(5) exp (Xt)(i(z), iv'(z)), with (v, iv) bounded and with Re A > 0, then (vj(z), wc(z)) is stable (3-5). On the other hand, if there is such a solution to Eq. 4, then (vc(z), wc(z)) is unstable and will change at an exponential rate of Re X in the direction of the unstable mode (v(z), iw'(z)). Note that Eq. 5 with A = 0 and (v, v) = (vs, w') is always a solution of Eq. 4, again because translates of (v¢, wc) are also solutions to Eq. 2. The time derivatives of Eq. 5 simplify, so for any A the function (v, wv) must satisfy

(v(z, t), w(z, t))

=

- cv' = (1 + A - b(z)/v4(0) + b(z - Z1)/v(z1))i + wc

= bv - (d + X)*,

(6)

a linear equation for all z with jump conditions on v' (z) across 0 and z1 . The search for values of A (eigenvalues) which support bounded solutions (eigenvectors) of Eq. 6 is very similar to the search for values of c which support impulse solutions to Eq. 2. If Re A > 0, the structure of the exponential functions which compose solutions again uniquely determines a solution (up to multiplication by a constant) bounded as z -0 -00. In particular, (v1, wv) on the interval (-o ,0) must simply be a multiple of the sole exponential factor (v+, w+) whose exponent has a positive real part. If such a solution (v, *) is an eigenvector, it must also be bounded as z +o, and so the coefficient of (v+ , w+) on the interval (z,, oo) (i.e., after the last jump) must be zero. We call that co-

efficient D(A). The function D(A) was originally defined and studied in a more general nerve equation setting by Evans (6). We use it to locate the eigenvalues X for Eq. 6 and we will need the following of its properties: (a) D(A) is complex analytic for all A, with Re A > 0. (b) D(A) = 0 if and only if A . (d) D(A) = D(X). isaneigenvaluefor Eq. 6. (c) ReD(A)-- oo as A 00, the numBecause D(A) is analytic with the above asymptotic behavior as A

ber of zeros with Re A > 0 (counting multiplicities) is given by the winding number of the image of the imaginary axis. The computational procedure used was as follows: select values for a, b, and d; find the values of c which support a solitary impulse solution; for each impulse compute D(A) for imaginary A to find the number of unstable modes; locate each eigenvalue X; compute each unstable eigenvector (IF, w). In the figures provided, b 0.2 and d 0.05, while the value of a was allowed to vary. Fig. 3 shows for a = 0.2 the image of the imaginary axis under D for both the

106

BIOPHYSICAL JOURNAL VOLUME 21 1978

fast cf, and slow, ca, velocities. The image for c, (Fig. 3 S) winds once, indicating the presence of one unstable eigenvector, while the image for cf (Fig. 3 F) does not wind, has no unstable eigenvector, and shows that (va, w,) for c = cf is a stable solution to Eq. 1. Fig. 4 shows in detail near the origin the image of the positive imaginary axis under D for several values of a and for both the slow (solid lines) and fast (dotted lines) velocities. Since D(X) = D(X), the image for corresponding A values is found by reflecting the given curve across the real axis. Again each slow impulse is seen to be unstable while each fast impulse is seen to be stable. The slow and fast curves will converge to a curve with a horizontal tangent at A = 0 as the increasing value a converges to the value which supports a unique impulse with velocity c,. That impulse has A = 0 as an eigenvalue of algebraic multiplicity two and so is called neutrally stable. Fig. 5 shows the graph of D(A) with (real) A > 0 for a = 0.24 and both c = c3, c = cf. The eigenvalue A > 0 is identified as the value on the slow (solid) graph where D(A) = 0. There is (necessarily) no positive eigenvalue for cf (dotted graph). A more

F

0.2602 0.275 tt 0.275

S

0.26 -0.5

0.5-0.5

03

~~~~~~~~0.21

0.5

FIGURE 3

FIGURE 4

FIGURE 3 The images in the complex plane of the imaginary axis as mapped by the function D(A) for the fast (F) and slow (S) solitary impulses. The loop in S indicates the presence of one unstable growth mode for the slow impulse. The absence of a loop in F shows that the fast impulse is stable. Parameter values used here are a = 0.2, b = 0.2, d = 0.05, c, = 0.4158, cfj

1.420. FIGURE 4 A detail of the complex plane near the origin showing the images of the positive imaginary axis under mapping by D(A) for various values of a and for both c, (solid lines) and Cf (dotted lines). Upon reflecting each curve across the real axis, the slow impulse curves show the characteristic loop seen in Fig. 3 S, while the fast impulse curves never have such a loop. The fixed parameter values here are b = 0.2, d = 0.05. The image for a = 0.2 and for c = cf is very close to that for a = 0.24 and is therefore not displayed.

J. A. FEROE Solitary Impulse Solutions of a Nerve Equation

107

0.3 4-

0.1

,-

~~

..'

~~~~~~~~~~~~~~~0 15

~~~~~~~~~~~~~~Az \/A

-,

X (Real)

-0.1

FIGURE 5

'..

FIGURE 6

FIGURE 5 The graphs of D(A) for positive (real) A for both the slow (solid graph) and fast (dotted graph) solitary impulses. The growth rate for the unstable mode of the slow impulse is approximately A = 1.433, at which value D(A) = 0. The graph of D(A) has no positive zero for the fast stable impulse. The parameter values here are a = 0.24, b = 0.2, c = 0.05, c. 0.4720, and 1.0779. FIGURE 6 The graph (solid line) of the unstable mode v(z) for the slow impulse vc(z) (dotted line). The direction of growth is both to change the width and height of the impulse. Parameter values here are a = 0.24, b = 0.2, d = 0.05, c. = 0.4720, and A = 1.433.

complete figure indicating the dependence of the unstable eigenvalue on the values of a appears in ref. 15. In Fig. 6 the unstable mode eigenvector (solid line);(z) for a = 0.24, c = c. is shown superimposed on vc(z) (dotted line). As one would expect, the initial direction of change is to increase or decrease both the height and the width of the impulse. These computations are representative of those for all parameter values for which computations have been carried out, and support the following proposition: Every fast solitary impulse solution of Eq. 1 is stable. Every slow solitary impulse solution of Eq. 1 is unstable with a single mode of instability. The unique solitary impulse with c = C. is neutrally stable. APPENDIX We study Eq. 2 in the expanded first order system

{VC

\WC

I=

0b

°

-

C

0O b/c -d/C \0/

II(vc - a)( 0

(7)

°

which we write as

V'

108

=

AV

-

h(V),

BIOPHYSICAL JOURNAL VOLUME 21 1978

by making the appropriate identifications. The constant matrix A has eigenvalues a 1, a2, a3 with corresponding eivenvectors Y1, Y2, Y3 where a I > 0, Re a2 < 0, Re a3 < 0, and a3 = a2 if a2 is complex. We define the eigenvectors by Y, = '(a,, 1, b/(d + aic)), and let E(z) be the fundamental matrix for Eq. 7, which has exp(aiz)Y, as its Pth column. If v,(z*) < a at some initial position z., the solution to Eq. 7 with initial condition V(z*) is locally V(z) = E(z Z*)E-'(O)V(z.), as long as v,(z) < a. For v,(z*) > a the solution to Eq. 7 is locally V(z) = E(z - z*)E-'(O)T(z.) + U, as long as v,(z) > a, where U = A -"((l,O,O), T(z.) = V(z*) U. At the transition points 0 and zI for a pulse solution (vc, wc) where vc(O) = vc(zi) = a, the fact that vc(z) - 0 as z |- co requires that only the first component of E-I(O)V(zo) be nonzero (in fact it requires that V(zo) = aY, ), and that the first component of E'-(O)V(z1) must be zero. Eq. 6 in the expanded first-order form is

{v = I1 0

0 b/c

0 -(d + A)/

v iw

(8)

for all z other than 0 and z, and subject to the jump conditions

O+ V

O= - v(0)/Vc(0) L'

+z =

V(Z1)/V'(Z1).

We write Eq. 8 as V' = A, V, by making the appropriate identification. If Re A > 0, the matrix A , will have eigenvalues f3, ,B2, 3 with corresponding eigenvectors X 1, X2, X3 where Re I > 0, Re B2 < 0, Rej33 < 0, and where we choose X, = '(fli, 1,b/(d + lic + A)), i = 1,2,3. Let M(z) be the fundamental matrix for Eq. 8 whose ith column is given by exp(,iz)X,. Solutions of Eq. 8 with initial condition V(z*) at z. are of the form V(z) = M(z - z.) M- (0)V(z), subject of course to the jump conditions across 0 and z1 . The condition that V be bounded demands that only the first component of M- (O)V(O -) be nonzero. Define B(A, z) to be the unique solution to Eq. 8 with B(A, 0 -) = X I. The unstable growth modes are precisely those A with Re A > 0 for which B(A,z) is bounded as z - +o, that is, those for which M'-(0)B(A, zt ) has a zero first component. The function D(A) is defined as the first component of

M'-(O)B(Aztz). This work was supported in part by National Science Foundation grant MCS77-00839. Receivedjor publication 19 August 1977.

REFERENCES 1. CASTEN, R., H. COHEN, and P. LAGERSTROM. 1975. Perturbation analysis of an approximation to Hodgkin-Huxley theory. J. Appl. Math. 32:365. 2. COOLEY, J. W., and F. A. DODGE. 1966. Digital computer solutions for excitation and propagation of the nerve impulse. Biophys. J. 6:583. 3. EVANS, J. W. 1972. Nerve axon equations. I. Linear approximations. Indiana Univ. Math. J. 21:877. 4. EVANS, J. W. 1972. Nerve axon equations. II. Stability at rest. Indiana Univ. Math. J. 22:75. 5. EVANS, J. W. 1972. Nerve axon equations. III. Stability of the nerve impulse. Indiana Univ. Math. J. 22:577. 6. EVANS, J. W. 1975. Nerve axon equations. IV. The stable and the unstable impulse. Indiana Univ. Math. J. 24:1169.

J. A. FEROE Solitary Impulse Solutions of a Nerve Equation

109

7. EVANS, J. W., and J. A. FEROE. 1978. Local stability theory of the nerve impulse. Math. Biosci. In press. 8. FITZHUGH, R. 1961. Impulses and physiological states in models of nerve membrane. Biophys. J. 1: 445. 9. FITZHUGH, R. 1969. Mathematical models of excitation and propagation in nerve. In Biological Engineering. H. P. Schwan, editor. McGraw-Hill Book Company, Inc., New York. 1-85. 10. HODGKIN, A. L., and A. F. HUXLEY. 1952. A quantitative description of membrane current and its application to condution and excitation in nerve. J. Physiol. (Lond.). 117:500. 11. HUXLEY, A. F. 1959. Can a nerve propagate a subthreshold disturbance? J. Physiol. (Lond.). 148:80P. 12. McKEAN, H. P. 1971. Nagumo's equation. Adv. Math. 4:209. 13. NAGUMO, J., S. ARIMoTo, and S. YOSHIZAWA. 1962. An active pulse transmission line simulating nerve axon. Proc. Inst. Radio Eng. 50:2061. 14. RINZEL, J. Spatial stability of traveling wave solutions of a nerve condution equation. Biophys. J. 15: 975. 15. RINZEL, J., and J. B. KELLER. 1973. Traveling wave solutions of a nerve condution equation. Biophys. J. 13:1313.

110

BIOPHYSICAL JOURNAL VOLUME 21 1978