Waves in excitable media

(1)

T.M.J. van Zalen

Waves in excitable media

Master thesis, defended on June 13, 2013 Thesis advisor: dr. S.C. Hille

Specialisation: Applied Mathematics

Mathematisch Instituut, Universiteit Leiden

(2)

(3)

Abstract

In this thesis we will consider waves in excitable media. These waves can be found in dierent biological settings, from signalling in the amoeba Dictyostelium discoideum to cardiac brillation and spiral galaxies. De- pending on the number of spatial dimensions, the waves can be of dierent types.

In one spatial dimension travelling waves and travelling pulses occur.

The existence and uniqueness of such waves will be proven. It follows that there there is a unique velocity, depending on the level of a so-called controller species, for which there exists a corresponding travelling wave.

In two spatial dimensions the spiral waves will be considered. It follows that the normal velocity of a spiral wave front is not only determined by the level of controller species in the front, but also by the curvature of the front.

(4)

1 Introduction on waves in excitable media

Excitable media have special characteristics: there is a stable rest state and small perturbations from the rest state are rapidly damped out. If the perturbation crosses a certain threshold, there will be an abrupt and substantial response. The medium will be refractory to further stimulation for a while until it recovers full excitability.

Waves in excitable media are found in for example the Belousov-Zhabotinskii reaction [4] and the cAMP signalling in the cellular slime mold Dictyostelium discoideum, which will be described in Chapter 2. Waves also occur on larger scale, such as the waves of infectious diseases that travel through populations, and a spiral galaxy, whose rotating arms can be treated as travelling waves of star formation in an excitable medium of interstellar gas and dust [9].

In biology it is not always a good thing for a system that relies on a faithful propagation of a signal to be taken over by a self-sustained pattern: spirals on the heart lead to cardiac brillation [8], spirals in the cortex may lead to epileptic seizures and spirals on the retina or visual cortex may cause hallucinations.

In excitable media there are dierent types of waves depending on the number of spatial dimensions. In one spatial dimension, we can distinguish between wave fronts and wave backs and between travelling pulses, which contains one wave front and wave back, and wave trains, which are sequences of wave fronts and wave backs and can be periodic. In Chapter 4 the travelling waves and travelling pulses in one dimension will be analysed in detail.

In two dimensions there are two characteristic patterns: the target patterns, which are created if the medium is suciently large so that more than one wave front can exist at the same time, and the rotating spiral waves. Target patterns are expanding circular waves and require a periodic source and cannot exist in a homogeneous non-oscillatory medium, while spiral waves do not require a periodic source for their existence, since they are typically self-sustained.

In Chapter 5 the spiral waves will be analysed in more detail.

Experimental studies showed that spiral waves are generated by spiral breakup, which was found to occur in models where the recovery variable diuses at a high rate. In [5] is was shown that the spiral breakup also can be caused by lateral instability of the wave front. A spiral wave can break up into a large number of small spirals. As showed in [14], in Dictyostelium these small spirals can spontaneously form multi-armed spirals in low excitable media. These multi- armed spirals are formed due to attraction of single spirals that rotate in the same direction.

In three dimensions the target patterns and rotating spiral waves generalize to expanding spherical waves and rotating scroll waves respectively. In this thesis the three dimensional cases will not be analysed. The interested reader may consult [12] and the references found there.

In this master thesis the focus is on understanding all details involved in the construction of travelling waves and pulses in one spatial dimension and spiral waves in two dimensions, as presented in a series of papers, [2], [3], [4] and [12].

We will not discuss the complete dynamics of the system of partial dierential equations (PDEs) in innite dimensional state space, such as stability of the waves, pulses and spirals. The construction of the latter 'reduces' the arguments to a setting of ordinary dierential equations (ODEs).

(6)

2 Dictyostelium discoideum: cAMP signalling

2.1 The biochemistry of the Martiel-Goldbeter model

In this chapter we will analyse the cyclic Adenosine 3'-5'-Monophosphate (cAMP) signalling system in Dictyostelium discoideum. The biological background of this cAMP signalling is the following, [13].

If the Dictyostelium discoideum amoebae are left to starve on an agar surface, they start to signal to each other by secreting cAMP. The individual cells receive this cAMP by binding extracellular cAMP to a membrane receptor. This stimulates the synthesis of cAMP from ATP by adenylate cyclase within the cell.

The newly synthesized cAMP is transported to the extracellular medium so that the chemical signal is amplied. However, the amplication is limited since in prolonged exposure to cAMP the membrane receptors become desensitized so that it no longer stimulates adenylate cyclase activity. In the absence of cAMP synthesis, the concentration of cAMP decreases by the action of extracellular phosphodiesterase, which hydrolyses cAMP.

Pulses of extracellular cAMP, secreted by the signalling amoebae, travel across the eld in the form of either expanding concentric circular waves or rotating spiral waves. As the waves pass periodically through the eld of inde- pendent amoebal cells, they stimulate a chemotactic movement of the amoebae towards the center of the pattern, which is the origin of the circles or the pivot point of the spiral. Eventually all the amoebae within the domain of a single pattern aggregate at the center to form a multicellular slug, which goes on to form a fruiting body.

The model for cAMP signalling (see Figure 2.1) in Dictyostelium corresponds to the sequence of reaction steps (a)-(e) found in (2.1), [6]. Although this is a hypothetical mechanism, there seems no other mechanism proposed for the cAMP signalling.

Figure 2.1: (a) Model for the cAMP signalling system of Dictyostelium discoideum. (b) The receptor box with the interconversions between the R and D state in the presence of P .

(7)

(a) R ^k_k¹₋₁ D (b) R + P âd₁¹RP (c) D + P â_d²₂ DP (d) RP ^k_k²₋₂ DP (e) 2RP + C â_d³₃ E

(f ) E + S ^da⁴4 ES →^k⁴E + Pi

(g) C + S ^da⁵₅CS →^k⁵ C + Pi

(h) Pi→^kⁱ (i) Pi→^k^t P →^k^e

(j) →^vⁱS →^k⁰ (2.1)

In reaction steps (2.1), R and D denote the active and desensitized form of the receptor respectively. Reactions (a)−(d) show the interconversions between the R and D states in the presence of extracellular cAMP (P ). In reaction (e), two molecules of active receptor-cAMP complex (RP ) bind to less-active adenylate cyclase (C) to form an enzyme complex of higher order activity (E).

Reactions (f) and (g) describe the synthesis of intracellular cAMP (Pi) from ATP (S) by both forms of adenylate cyclase. In reaction (h) part of the Pi is hydrolysed by an intracellular phosphodiesterase, and in reaction (i) part of Pi

is transported across the plasma membrane into the extracellular medium where it is hydrolysed by the membranebound and extracellular forms of this enzyme.

Finally, it is assumed that the substrate ATP is synthesized at a constant rate and used in a number of reactions besides that catalysed by adenylate cyclase.

2.2 The mathematical models

Considering the time evolution of the areal densities of the species, we are able to nd a system of dierential equations corresponding to the reactions in (2.1).

First of all, we need to introduce some notations:

Denition 1 (Areal density). The areal density of a membrane bound molecular species is given by [·]A.

If Am is the area of the membrane, Ve the extracellular volume and Vi the intracellular volume, we can also dene, following [6]

Denition 2 (Concentration). The concentration of a membrane bound intracellular species is given by

[·_i] := [·_i]_AA_m Vi

, (2.2)

and the concentration of a membrane bound extracellular species is given by

[·_e] := [·_e]_AA_m Ve

. (2.3)

(8)

Here, the membrane bound intracellular species are C, E, CS, ES, S and P_iand the extracellular species are R, P , D, RP and DP (see Figure 2.1). For simplicity we will omit the subscripts i and e in the rest of the chapter.

Further, we will write the reaction rates for the dierential equations of the areal densities as k and d for the simple rates, and ˆa and ˆv for the other rates.

If we look at the time evolution of the areal density of the species R, we nd the following dierential equation:

d[R]_A

dt = −k₁[D]_A+ k₋₁[R]_A− ˆa₁[R]_A[P ]_A+ d₁[RP ]_A. (2.4) We will make use of two conservation relations for the receptor and for adenylate cyclase:

RT = [R]A+ [D]A+ [RP ]A+ [DP ]A+ (2/h)([E]A+ [ES]A) (2.5)

E_T = [E]_A+ [ES]_A+ [C]_A+ [CS]_A, (2.6)

with h a dilution factor.

If we now substitute ρ = [R]/RT, δ = [D]/RT, γ = [P ]/RT and x = [RP ]/RT

in (2.4), we nd the following dierential equation:

dρ

dt = −k1δ + L1k1ρ − d1

ˆ a1

a1

Ve

Am

ργ + d1x, (2.7)

and if we take

ˆ a1= a1

Am

Ve

, we nd the equation

dρ

dt = k1(−δ + L1ρ) + d1(−ργ + x). (2.8) The denition and the values used for the calculations of the parameters can be found in Table 2.1, used in [11] for numerical simulations. Note that the parameters are dimensionless, except the time-scale, which is in minutes, and the space scale, which is in mm.

The reason that [6] uses concentrations rather than surface densities for the membrane bound molecular species, may be that the authors seem to link experimental data concerning total concentration of these compounds in a reaction vessel.

In this way we can derive the dierential equations for the time evolution of the concentration of the other species, which can be found in Appendix A in [6].

However, we had to make some adjustments and supplements to this article in order to derive the dierential equations:

c = [C]/ET, cs = [CS]/ET, y = [DP ]/RT

β = [Pi]/KR, γ = [P ]/KR, K_m⁰ = (d5+ k5)/a5.

If we use the method described above to derive the system of equations, we also

nd the reaction rates ˆ

a2= a2

A_m

, ˆa3= a3

A²_m

, ˆa4= a4

A_m

, ˆa5= a5

A_m

, ˆvi= vi

V_i .

(9)

The system of dierential equations given in Appendix A of [6] can be reduced to a four-variable system, under the assumptions that the dierential equations contain both fast binding and slow modication terms, so the following inequality on the rate constants holds

(k1, k₋₁, k2, k₋₂, ki, kt, ke, σ, k⁰) (a1, d1, a2, d2, a3, d3, a4, d4, a5, d5), and that the dierential equations for the fastest variables (the total fraction of the receptor forms bound to cAMP (x + y), c, e and es) reduce to algebraic equations corresponding to the quasi-steady-state hypothesis for these receptor and enzyme forms. The evolution equations for the remaining slow variables can now be transformed, yielding the four-variable system of dierential equations (2.9).

dρT

dt = −f₁(γ)ρ_T + f₂(γ)(1 − ρ_T) dα

dt = v − k⁰α − σΦ(ρT, γ, α) dβ

dt = qσΦ(ρT, γ, α) − (ki+ kt)β dγ

dt = (k_tβ/h) − k_eγ (2.9)

with

f₁(γ) = k₁+ k₂γ

1 + γ , f₂(γ) = k₁L₁+ k₂L₂cγ 1 + cγ Φ(ρ_T, γ, α) = α(λθ + εY²)

1 + αθ + εY²(1 + α), Y = ρTγ 1 + γ.

Here, ρT denotes the total fraction of receptor in the active state, β and γ denote the intracellular and extracellular concentrations of cAMP divided by the dissociation constant KR, and α is the intracellular ATP concentration divided by the Michaelis constant Kmof reaction (f), so that Km= ^d⁴_a^+k⁴

4 . Experiments have indicated that the level of ATP can be kept constant, since the absence of signicant variation in ATP can be accounted for by the model. Hence we can consider α as a parameter, and we are able to further reduce the number of variables. The dynamics of the cAMP signalling system is then governed by the three dierential equations

dρ_T

dt = −f₁(γ)ρ_T + f₂(γ)(1 − ρ_T) dβ

dt = qσΦ(ρT, γ, α) − (ki+ kt)β dγ

dt = (ktβ/h) − keγ (2.10)

We can make time dimensionless by taking t = k1· t, since we consider the time scale in which the receptor will desensitize. Substituting this into

(10)

equations (2.10) results in the dierential equations dρT

dt = −f₁(γ)ρT+ f₂(γ)(1 − ρT) (2.11) ε⁰dβ

dt = s1Φ(ρT, γ) − β (2.12)

ε⁰⁰dγ

dt = s2β − γ (2.13)

where

f₁(γ) = 1 + κγ

1 + γ , f₂=L1+ κL2cγ 1 + cγ Φ(ρT, γ) = λ1+ Y²

λ₂+ Y².

Name Denition Set A Set B Set C Set D Set E

L₁ k₋₁/k₁ 10 = = = =

L₂ k₋₂/k₂ 0.005 0.005 0.005 0.005 0.0005

κ k₂/k₁ 18.5 = = = =

c K_R/K_D 10 10 10 10 45

α [ATP]/Km 3 = = = =

λ1 λθ

ε 10⁻⁴ 10⁻³ 10⁻³ 10⁻³ 6.7 · 10⁻⁴

λ₂ _(1+α)ε^1+αθ 0.26 2.4 2.4 2.4 1.0

s₁ _k^qσ

i+kt

α

1+α 690 950 950 360 80

s2 kt/keh 0.033 0.05 0.05 0.13 0.35

s s1s2 23 47 47 47 28

ε⁰ k1/(ki+ kt) 0.014 0.019 0.019 0.005 0.01

ε⁰⁰ k₁/k_e 0.0067 0.01 0.01 0.01 0.024

Time-scale 1/k₁ 28 28 8.3 28 17

Space-scale (keD)¹²/k1 10 8.2 4.5 8.2 4.1 Table 2.1: Model parameters as occuring in [6] (A-C) and [11] (D-E).

When all four parameter sets assume the same parameter value, the symbol

0 =⁰ is used. Note that set A is used in [6] to model cAMP oscillations in a well-stirred cell suspension. Set B is used in [6] to model cAMP signal-relaying in a well-stirred cell suspension. Set C is used in the paper by Tyson et al. to calculate spiral waves in the full three-component model. Sets D and E are used in [11] to calculate spiral waves in the two-component model.

Equations (2.11)-(2.13) describe a homogeneous distribution of the species.

However, if we allow the species to diuse freely, we must include diusion. The equations describing the membrane-bound receptor and the intracellular cAMP are unchanged since these species are immobile, since located within cells. The extracellular cAMP can freely diuse, so it must be described by a dimensionless equation in the form of a reaction-diusion equation (for simplicity we will write t for the scaled time parameter)

∂γ

∂t = ε⁰⁰∆γ + 1

ε⁰⁰F (β, γ),

(11)

with F (β, γ) the right-hand-side of equation (2.13), and ∆ the two-dimensional Laplacian. The spatial variables are scaled so that the diusion coecient is equal to ε⁰⁰, that is

(x, y) = kt

√k_eD · (x, y). (2.14)

Here x, y are the dimensional spatial variables and D is the diusion coecient of cAMP. The length scale used in (2.14) is the length scale of exponential decay of the steady-state solution of diusive species on (0, ∞), with diusivity D, that degrade with rate ke and enter at 0 with rate kt. According to [11], we choose D = 0.024mm²/min. The biological reason to take this space scale is that we consider the combination of diusion D and degradation ke of the extracellular cAMP.

Analysis of this Martiel-Goldbeter model would be simplied if we could eliminate the equation for β. Such simplication would be possible if ε⁰ ε⁰⁰, so that a quasi steady-state of β can occur. This quasi-steady state will be β = s1Φ(ρT, γ). In that case we can do a phase plane analysis of the relay and oscillations. For the parameter sets in [11] parameter sets D and E can be used, so that it is justied to take ε⁰→ 0with ε⁰⁰ nite. In that case the model equations become

∂γ

∂t = ε⁰⁰∆γ + 1

ε⁰⁰[sΦ(ρ_T, γ) − γ]

∂ρT

∂t = −f₁(γ)ρT + f₂(γ)(1 − ρT), (2.15) where s = s1s₂.

Considering parameter sets D and E, we can compute the nullclines of the ODE part of system (2.15). The γ-nullcline is

0 = sΦ(ρT, γ) − γ

0 = s

λ1+_ρ

Tγ 1+γ

2

λ₂+

ρTγ 1+γ

2 − γ

s λ1+ ρTγ 1 + γ

²!

= γ λ2+ ρTγ 1 + γ

²!

ρT(γ) = ±1 + γ γ

s

λ2γ − sλ1

s − γ . (2.16)

Since we only consider positive values of ρT, the only nullcline we will consider is the positive solution branch.

From equation (2.16) we can conclude that the γ-nullcline does not have points with γ < ^sλ_λ₂¹ or γ > s.

The ρT-nullcline is

0 = −f₁(γ)ρT+ f₂(γ)(1 − ρT) ρT(γ) = f₂(γ)

f₁(γ) + f₂(γ)

= (L₁+ κL₂cγ)(1 + γ)

. (2.17)

(12)

We can also compute the maxima of the γ-nullcline (2.16). First we have to determine ^∂ρ_∂γ^T

∂ρT

∂γ =

"

1 + γ γ

s

λ2γ − sλ1

s − γ

#⁰

(2.18)

= (γ + γ²)(λ2− λ1)s − 2(λ2γ − sλ1)(s − γ) 2γ²(s − γ)²q

λ2γ−sλ1

s−γ

(2.19)

= γ²(2λ2+ s(λ2− λ1)) − γs(λ2+ 3λ1) + 2s²λ1

2γ²(s − γ)²

qλ₂γ−sλ₁ s−γ

(2.20)

Solving ^∂ρ_∂γ^T = 0for γ results in

0 = γ²(2λ2+ s(λ2− λ1) − γs(λ2+ 3λ1) + 2s²λ1 (2.21) γ_± = s[3λ1+ λ2±p(λ2− λ1)(λ2− (8s + 9)λ1]

2s(λ₂− λ₁) + 4λ₂ (2.22)

Equation (2.22) is real if λ2 < λ1 or λ2> (8s + 9)λ1. According to the model parameters for λ1 and λ2, see Table 2.1, λ1 λ2, hence the only possibility for the maxima of the nullcline to exists, is that λ2> (8s + 9)λ₁. In that case, the sign of the denominator is s(λ2− λ₁) + 2λ₂> 0. Hence, we can conclude that for γ < γ−, ρT is increasing, for γ−< γ < γ₊, ρT is decreasing and for γ > γ+, ρ_T is increasing again, so that this nullcline is N-shaped, which can be seen in

gure 2.2.

The nullclines of system (2.15), with parameter set D, are illustrated in Figure 3.1. Note that in this gure γ is plotted against 1−ρT, since the recovery variable in this case is the re-sensitization of the receptor ρT.

Figure 2.2: Nullclines of system (2.15) for parameter set D.

(13)

3 A reduced two-state model

The pair of dierential equations given in system (2.15), is of the form

ε∂u

∂t = ε²∆u + f (u, v),

∂v

∂t = g(u, v), (3.1)

with

u = γ, v = 1−ρ_T, f (u, v) = sΦ(1−ρ_T, γ)−γ, g(u, v) = −f₁(γ)(1−ρ_T)+f₂(γ)ρ_T. In the rest of this thesis we will consider a pair of reaction-diusion equations of the form

ε∂u

∂t = ε²∆u + f (u, v),

∂v

∂t = εδ∆v + g(u, v),

(3.2)

with u = u(x1, . . . , xn, t) and v = v(x1, . . . , xn, t), 0 < ε 1 the ratio of the rates of reaction of the quantities u and v, ∆ = Pⁿi=1

∂

∂²x_i the Laplacian operator with n the dimensions, and δ = _D^D^v_u the ratio of diusion coecients of the two components u and v. u can be seen as the trigger or excitation variable and v as the recovery variable.

We will consider in particular two cases of system (3.2): δ = 0 and δ = 1.

If δ = 0 the recovery variable v does not spread in space. This is characteristic in activity waves in neuromuscular tissue where v represents the local permeability of a membrane to transmembrane ionic currents. The two-dimensional case corresponding to δ = 0 has also been considered in [3]. Moreover, Dictyostelium corresponds to this case.

If δ = 1 both components u and v diuse equally well through space. This is characteristic for the Belousov-Zhabotinskii reaction. The two-dimensional case corresponding to δ = 1 has also been considered in [4].

The functions f(u, v) and g(u, v) describe the non-linear kinetics of system (3.2). We make the following assumptions on the function f(u, v), following [2]:

(AF1) f is C¹ and the nullcline f(u, v) = 0 is cubic-shaped, hence there are three solution branches to f(u(v), v) = 0. Denote them by u = h±(v)and u = h0(v), with h−(v) ≤ h0(v) ≤ h+(v). Here h− is dened on (v−, ∞), h0 on [v−, v+]and h+ on (−∞, v+). We require that v−< v+.

(AF2) ^∂f_∂u(h_±(v), v) < 0, so these solution branches are stable, and^∂f_∂u(h0(v), v) ≥ 0. Therefore this solution branch is unstable.

(AF3) Above the line f(u, v) = 0, f(u, v) < 0 and below this line f(u, v) > 0.

(AF4) ^∂f_∂v(u, v) ≤ 0on u ∈ [h−(v), h+(v)]for each v ∈ [v−, v+].

Example 1. 1. Let f(u, v) = P3(u)−v, with P3(u)a third order polynomial of u. Then ^∂f_∂v = −1 < 0for all v.

(14)

2. For Dictyostelium discoideum we found the following function f(u, v):

f (u, v) = sλ1(1 + u)²+ (1 − v)²u² λ2(1 + u)²+ (1 − v)²u² − u.

If we calculate ^∂f_∂v we obtain:

∂f

∂v(u, v) = s−2(1 − v)u²(1 + u)²(λ2− λ1) (λ₂(1 + u)²+ (1 − v)²u²)² . Since λ2> λ1, u > 0 and 1 − v > 0, it yields that

∂f

∂v(u, v) < 0.

Because of the Implicit Function Theorem, h± are C¹functions. Dene

If(v) :=

Z h₊(v) h−(v)

f (z, v) dz. (3.3)

Proposition 1. Assume (AF1)(AF3). If(v)has the following properties:

(i) If ∈ C¹([v₋, v+]).

(ii) There exists v^∗±: v−< v₋^∗ ≤ v₊^∗ < v+ such that If(v0) > 0 for v−≤ v0<

v^∗₋ and If(v0) < 0 for v+^∗ < v0≤ v+.

(iii) If in addition (AF4) holds, then If⁰(v) < 0on (v−, v₊). Consequently there exists a unique v^∗∈ (v−, v₊) such that If(v^∗) = 0and v^∗−= v^∗₊= v^∗. Proof. (i) Because f is C¹, h± are C¹functions. The Fundamental Theorem

of Calculus then yields that If is C¹. (ii) This follows directly from (AF3).

(iii) Because of the denition of h±, one has If⁰(v) =Rh₊(v) h−(v)

∂f

∂v(z, v) dz. Since (AF4) holds, ^∂f_∂v(u, v) < 0and thus If⁰(v) < 0. From (i) and (iii) follows directly that there should be a v^∗∈ (v₋, v+)such that If(v^∗) = 0.

Following [2], the assumptions made on g(u, v) are:

(AG1) The nullcline g(u, v) = 0 is monotone increasing.

(AG2) g(u, v) = 0 has precisely one intersection with the curve f(u, v) = 0. This intersection is the unique steady state solution and is located at the branch u = h₋(v). It is denoted by (us, v_s), such that v−< v_s< v^∗.

(AG3) g(u, v) < 0 above (to the left) the line g(u, v) = 0 and g(u, v) > 0 below (to the right) this line.

(15)

The characteristic setting of these nullclines is illustrated in Figure 3.1.

f (u, v) = 0 g(u, v) = 0

u v

f < 0 g < 0

f > 0 g > 0 h−(v)

h0(v) h+(v)

us vs v+

v−

Figure 3.1: Typical nullclines f(u, v) = 0 and g(u, v) = 0 for the excitable system (3.2).

The structures for f and g as described above are typical for excitable sys- tems. A system is excitable if a perturbation over a certain threshold triggers a large excursion from the steady state before it eventually returns to the steady state. In other words, a system is excitable if, with a sucient initial stimulus, a pulse can be initiated which propagates throughout the medium.

The excursion can be described as follows. The excitation variable u increases rapidly until the trajectory approaches the branch f(h+(v), v). At this point the trajectory moves slowly up the nullcline as the recovery variable v increases.

When the trajectory reaches the top of f(h+(v), v), the trajectory moves to the the branch f(h−(v), v) as u rapidly decreases, which is followed by a slow decrease of v back to the rest state (us, vs).

Although a pulse may be initiated, it may slow down, reverse direction, and eventually annihilate itself. In that case it will not propagate throughout the medium.

The form of f and g is such that the system of reaction-diusion equations (3.2) has a so-called invariant region (cf. [10]), in particular Theorem 4.3 and Example 6.C. That is, there exists a (necessarily) rectangular region R in R²+, containing 0, for which the vector eld associated to the reaction part in (3.2) points inward everywhere. Then, a solution (u(t), v(t)) to (3.2) with initial condition (u0, v0) such that (u0(x), v0(x)) ∈ R for all x will have (u(t), v(t)) ∈ Rfor all t ≥ 0. R contains the steady state (us, vs). Thus bound- edness of solutions is assumed.

(16)

4 Waves in one dimension

Throughout this chapter we assume that the properties of f and g, respectively (AF1)(AF4) and (AG1)(AG3), hold.

4.1 Travelling waves

First we will consider (3.2) in one spatial dimension, so that ∆ = _∂x^∂²2. This system can be treated using singular perturbation theory. Since the system has a small parameter ε, we can consider two regions, which contains a wave front.

First of all, in the slowly varying region, when ε = 0, we nd

f (u, v) = 0, ∂v

∂t = g(u, v).

Since f(u, v) = 0 has two branches of solutions that are stable and one solution unstable to perturbations o the branch, the outer solution to (3.2) is

∂v

∂t = g(h±(v), v). (4.1)

We are looking for travelling waves so we have to nd solutions to system (3.2) of xed speed. The solutions will be of the form

u(x, t) = u(z), v(x, t) = v(z), z = x − ct, with c the wave speed dependent on v.

In that case we can write system (3.2) as

ε²u_zz+ εcu_z+ f (u, v) = 0

εδvzz+ cvz+ g(u, v) = 0. (4.2) The O(1) equations, which is the lowest order in ε obtained by setting ε = 0 in (4.2), are

f (u, v) = 0, cv_z+ g(u, v) = 0.

4.1.1 Transition from h− to h+

Consider the transition from h−(v) to h+(v) (the so-called "up-jump"). If we suppose that for z < 0, u = h+(v)and for z > 0, u = h−(v), the solution u(z) is no longer continuous at z = 0. In that case diusion becomes important, i.e.

the terms ε²uzz and εcuz, and the regions of the solution branches are patched together by a boundary layer of width O(ε).

Therefore we introduce a boundary layer coordinate (or a stretching transfor- mation)

ξ = z ε.

When using this boundary layer coordinate, system (4.2) becomes u⁰⁰+ cu⁰+ f (u, v) = 0

δv⁰⁰+ cv⁰+ εg(u, v) = 0, (4.3)

(17)

with u⁰⁰= u_ξξ, u⁰= u_ξ etcetera.

The boundary layer should have the correct orientation, so limξ→−∞u(ξ) = h+(v)and limξ→∞u(ξ) = h₋(v), and c > 0, see Figure 4.1.

To lowest order in ε and requiring that v is bounded for all ξ, we have to solve the equations

u⁰⁰+ cu⁰+ f (u, v) = 0

δv⁰⁰+ cv⁰ = 0, (4.4)

with

lim

ξ→∓∞u(ξ) = h_±(v), lim

ξ→±∞u⁰(ξ) = 0.

Solving the second equation of (4.4), we nd the general solution v(ξ) = k1+ k2e⁻^c^δ^ξ,

and since v is bounded, see Chapter 3, for all ξ, we nd that k2 = 0 and v(ξ) ≡ k1, with k1 some constant. For simplicity we will write k1 = v0, with v0∈ [v−, v^∗].

Therefore, system (4.4) can be rewritten into

u⁰⁰+ cu⁰+ f (u, v) = 0, v = v₀, (4.5) with

lim

ξ→∓∞u(ξ) = h_±(v₀) lim

ξ→±∞u⁰(ξ) = 0. (4.6)

This transition is sketched in Figure 4.1.

h−(v0) h+(v0)

ξ u

|c|

Figure 4.1: Situation sketch for the transition from h−(v₀)to h+(v₀). System (4.5)- (4.6) is a non-linear eigenvalue problem for wave speed c parametrized by v0. We shall now discuss how c for which (4.5)- (4.6) has a solution, depends on v0.

Lemma 1. Let v−≤ v0≤ v+. Suppose c is such that (4.5)- (4.6) has a solution u(ξ), then

c = cu(v0) =

Rh+(v0)

h−(v0)f (u, v0) dz R∞

−∞(u(ξ)⁰)²dξ . (4.7) Proof. Consider

u⁰⁰+ cu⁰+ f (u, v₀) = 0, (4.8)

(18)

with limξ→∓∞u(ξ) = h_±(v₀)and limξ→±∞u⁰(ξ) = 0. Multiplying (4.8) by u⁰ and integrating from ξ = −∞ to ξ = ∞, we nd

Z ∞

−∞

u⁰⁰(ξ)u⁰(ξ) + c(u⁰(ξ))²+ f (u(ξ), v0)u⁰(ξ) dξ = 0 1

2(u⁰(ξ))²

∞

−∞

+ c Z ∞

−∞

(u⁰(ξ))²dξ + Z ∞

−∞

f (u(ξ), v0)u⁰(ξ) dξ = 0. (4.9) Since limξ→±∞u⁰(ξ) = 0,

1

2(u⁰(ξ))²

∞

−∞

= 0. (4.10)

Since u = u(ξ), we have du = u⁰dξ. Therefore, we can write Z ∞

−∞

f (u(ξ), v₀)u⁰(ξ) dξ =

Z u(∞) u(−∞)

f (z, v₀) dz

=

Z h₋(v₀) h+(v0)

f (z, v0) dz

= −

Z h₊(v₀) h−(v0)

f (z, v0) dz. (4.11) Substituting (4.10) and (4.11) into (4.9), we nd

c =

Rh+(v0)

h−(v₀) f (z, v0) dz R∞

−∞(u⁰(ξ))²dξ .

If a solution u to (4.5)- (4.6) exists for v− ≤ v₀≤ v₊, then c has the same sign as If(v0)(cf. equation (3.3). In particular

Corollary 1. If v−≤ v0< v^∗, then c > 0. If v^∗< v₀≤ v+, then c < 0.

Proof. From Proposition 1 we know that Z h₊(v₀)

h−(v₀)

f (u, v0) du > 0

for v− ≤ v0< v^∗. Combining this with (4.7) it follows that c > 0.

For v^∗< v₀≤ v₊ Proposition 1 states that Z h₊(v₀)

h−(v₀)

f (u, v0) < 0,

and thus c < 0. This corresponds with the down-jump in Section 4.1.2.

The properties of the wave speed c as a function of v0 is given in Figure 4.2.

To show that there exists a unique solution to (4.5)- (4.6), hence a unique c(v0) > 0 for all v− < v0< v^∗, we use the approach of [1], Chapter 4.2.1 (c).

Write u⁰ = w. Then (4.5) transforms into a rst-order system

u⁰ = w

w⁰= −cw − f (u, v ), (4.12)

(19)

v0 c

cmax

cmin

v− v∗ v+

Figure 4.2: Wave speed c against v0.

with

lim

ξ→∓∞(u, w) = (h_±(v0), 0).

Note that (u, w) = (h+(v₀), 0) and (u, w) = (h−(v₀), 0) are critical points for system (4.12). Thus we are looking for c such that there exists a heteroclinic orbit connecting (h+(v₀), 0) to h−(v₀), 0). First we determine the eigenvalues of the corresponding linearisation. The Jacobian matrix corresponding to the right hand side of system (4.12) is

J (u, w) =

0 1

−f⁰(u, v0) −c

. The linearisation around (h+(v₀), 0)then equals

u⁰ w⁰

= J (h+(v0), 0)

u − h₊(v₀) w

. Evaluating J(u, w) in the critical point (h+(v0), 0)gives

J (u, w) =

0 1

−f⁰(h+(v0), v0) −c

. (4.13)

The eigenvalues corresponding to (4.13) are

λ^±₁ =−c ±pc²− 4f⁰(h+(v0), v0)

2 . (4.14)

The eigenvalues for (h−(v0), 0)can be computed in a similar way, which results in

λ^±₀ =−c ±pc²− 4f⁰(h₋(v₀), v₀)

2 . (4.15)

Here the index 0 refers to the critical point x^∗0:= (h₋(v₀), 0)and 1 to the point x^∗₁:= (h₊(v₀), 0).

Since f⁰(u, v₀) < 0 for u = h±(v₀), λ^±0 and λ^±1 are real and both have dierent signs. The critical points are therefore saddle points. Hence, the unstable manifold W^u(x^∗₁) of x^∗1 corresponds to λ⁺1 and the stable manifold W^s(x^∗₀)of x^∗0 corresponds to λ⁻0. From this we can conclude that an unstable curve leaves x^∗ and a stable curve enters x^∗. Moreover W^u(x^∗)is tangent to

(20)

the line w = λ⁺1(u − h₊(v₀)) at x^∗1 and W^s(x^∗₀) is tangent to the line w = λ⁻₀(u − h₋(v₀))at x^∗0.

In order to nd a heteroclinic orbit connecting x^∗1and x^∗0, we have to nd c > 0 such that W^u(x^∗₁) = W^s(x^∗₀), for {u > 0, w < 0}.

Figure 4.3 shows a situation of the unstable manifold W^u(x^∗₁)and the stable manifold W^s(x^∗₀). Let Lεdenote the vertical line through (h0(v0) + ε, 0), with ε ≥ 0. For v− < v0 < v^∗ it is allowed to take ε = 0. In this case we write L instead of Lε. If v0 = v₋ then the critical point (h0(v₋), 0)coincides with x^∗0, and we must take ε > 0.

w

u

x∗0 (h0(v0), 0) x∗1

Lε W s (x∗0 )

W u (x∗1 ) w1(c)ˆ

w0(c)ˆ

Figure 4.3: Situation sketch for W^u and W^s.

Lemma 2. If c > 0, then there exist ˆw₀(c) < 0 and ˆw₁(c) < 0 such that W^u(x^∗₁) ∩ Lε= {(h0(v0) + ε, ˆw1(c))}, W^s(x^∗₀) ∩ Lε= {(h0(v0) + ε, ˆw0(c))}.

Proof. The unstable curve enters the region to the right of line Lεand cannot exit through the bottom, top or left hand side. Using (4.12) we deduce that W^u(x^∗₁)must exit this region through the line Lεat a point (h0(v₀) + ε, ˆw₁(c)). Similarly, we can argue that W^s(x^∗₀)must cross Lεat a point (h0(v₀)+ε, ˆw₀(c)). Moreover, the intersections consist of only one single point. Since u⁰ = w < 0 W^u(x^∗₁) cannot re-enter Lε and W^s(x^∗₀) cannot escape line Lε on the right side.

Lemma 3. If c = 0, then trajectories of (4.12) are contained in level sets of E(u, w) := ^w₂² +Ru

h−(v0)f (z, v₀) dz.

Proof. If E : R²→ R is dierentiable and is a constant of motion, then d

dtE(u(t), w(t)) = 0.

Hence

∂E

∂u du dt +∂E

∂w dw

dt = 0.

Consequently, since c = 0,

w∂E

∂u = f (u, v0)∂E

∂w. Thus we may take ^∂E_∂u = f (u, v0)and ^∂E_∂w = w. Then

E(u, w) =w² 2 +

Z u h−(v₀)

f (z, v₀) dz.

(21)

The level curves of E are of the form as illustrated in Figure 4.4.

u w

w1(0)ˆ w0(0)ˆ

Figure 4.4: Level curves of E(u, w).

At x^∗1 the level curves behave as shown in Figure 4.5.

x∗1

λ+ 1 λ−

1

Figure 4.5: Behaviour of the level curves at x^∗1.

Lemma 4. Let v−≤ v0< v^∗. If c = 0, it follows that ˆ

w₀(0) > ˆw₁(0). (4.16) Proof. If there is a heteroclinic orbit at c = 0, then x^∗0 and x^∗1 are on the same level set of E(u, v). In that case we have

E(h₋(v), 0) = E(h+(v), 0)

0 =

Z h+(v) h−(v)

f (z, v) dz. (4.17)

From equation (4.17) follows that if there exists a heteroclinic orbit with c = 0, then If(v) =Rh₊(v)

h−(v)f (z, v) dz = 0. From Proposition 1 follows that c(v) = 0 implies that v = v^∗.

(22)

On one hand we have

E(L, ˆw0(0)) = E(h−(v0), 0) = 0, (4.18) E(L, ˆw1(0)) = E(h+(v0), 0) =

Z h₊(v₀) h−(v0)

f (z, v0) dz, (4.19) while on the other hand

E(L, ˆw₁(0)) = wˆ₁(0)²

2 +

Z L h−(v₀)

f (z, v₀) dz, (4.20)

E(L, ˆw₀(0)) = wˆ₀(0)²

2 +

Z L h−(v₀)

f (z, v₀) dz. (4.21) Combining equations (4.18) and (4.21) gives

ˆ w₀(0)²

2 = −

Z L h−(v₀)

f (z, v0) dz. (4.22)

Note that − R_h^L₋_(v₀₎f (z, v₀) dz > 0. Combining equations (4.19) and (4.20) gives ˆ

w1(0)²

2 =

Z h₊(v₀) L

f (z, v0) dz, (4.23)

with R_L^h⁺^(v⁰⁾f (z, v0) dz > 0.

Since v−≤ v0< v^∗, If(v0) > 0, see Proposition 1. Thus ˆ

w₀(0)²

2 = −

Z L h−(v₀)

f (z, v₀) dz <

Z h+(v0) L

f (z, v₀) dz = wˆ₁(0)² 2 . It follows that | ˆw0(0)| < | ˆw1(0)|, and since ˆwi(0) < 0, ˆw0(0) > ˆw1(0). Lemma 5. Let v−≤ v0< v^∗. Then there exists cv₀ such that

ˆ

w0(c) < ˆw1(c), (4.24) provided c > cv₀.

Proof. We will show that ˆw₀(c) < ˆw₀(0) and ˆw₁(c) > ˆw₀(0) for c suciently large. Recall system (4.12).

The unstable manifold W^u(x^∗₁)is tangent to w = λ⁺1(c)(u − h+(v0)), with

λ⁺₁(c) =

−c +q

c²− 4^∂f_∂u(h₊(v₀), 0)

2 > 0.

For c large we have

c→∞lim λ⁺₁(c) = lim

c→∞

c²−q

c²− 4^∂f_∂u(h+(v0), v0) 2

1

−c −q

c²− 4^∂f_∂u(h+(v0), v0)

= lim

c→∞

2^∂f_∂u(h+(v0), v0)

−c − q

c²− 4^∂f_∂u(h₊(v₀), v₀)

= 0. (4.25)

(23)

Let β > 0. There exists cv0 > 0 such that |λ⁺1(c)| < β for all c ≥ cv0. If u ∈ (h₀(v₀), h₊(v₀)), f(u, v0) > 0. Along the line segment Tβ:= {h₀(v₀) < u <

h+(v0), w = β(u − h+(v0))}, see Figure 4.6, we have w⁰

u⁰ = −cw − f (u, v₀)

w = −c − f (u, v₀)

β(u − h+(v0))= −c + 1 β

f (u, v₀) h+(v0) − u

. Since

f (u,v₀) h+(v0)−u

≤ K, with K = supu∈[h₀(v₀),h₊(v₀)]

∂f

∂u(u, v₀)

, we see that w⁰

u⁰ ≤ −c +K β < |β|,

on Tβ, provided c > 0 is large enough. Therefore W^u(x^∗₁)cannot exit S through the u-axis or Tβ if we choose cβ is suciently positive and we conclude that

ˆ

w₁(c) > β(h₀(v₀) − h₊(v₀)), if c > cv0.

Hence for β > 0 large enough, such that |λ⁺1| < β, and c > 0 suciently positive, ˆ

w1(c) > ˆw0(0), (4.26) for all c > cβ.

w

u

(h0(v0), 0) x∗1

L Tβ S

λ+ 1

Figure 4.6: Situation sketch for region S.

The stable manifold W^s(x^∗₀)is tangent to w = λ⁻0(c)(u − h₋(v0)), with

λ⁻₀(c) =

−c − q

c²− 4^∂f_∂u(h−(v0), v0)

2 < 0,

with

c→∞lim λ⁻₀(c) = −∞, so ˆw₀(c)is unbounded for large c.

Now let β⁰ < 0. If u ∈ (h−(v₀), h₀(v₀))we have f(u, v0) < 0. Along the line segment Tβ⁰ := {h₋(v₀) < u < h₀(v₀) : w = β⁰(u − h₋(v₀))}, see Figure 4.7, we have

w⁰

u⁰ = −c −f (u, v₀) w

= −c +|f (u, v0)|

w

= −c + |f (u, v0)|

0

(24)

Since

f (u,v₀) u−h−(v0)

≤ M, with M = supu∈[h−(v₀),h₀(v₀)]

∂f

∂u(u, v0)

, we have w⁰

u⁰ ≤ −c +M β⁰ < β⁰,

on Tβ⁰, provided c > 0 is large enough. Therefore we nd that ˆ

w0(c) < β⁰(h0(v0) − h−(v0)),

if cβ⁰ is suciently large. Hence, W^s(x^∗₀)cannot exit S⁰ through Tβ⁰. Hence, for all β⁰< 0 there exists a cβ⁰ > 0such that

ˆ

w0(c) ≤ β⁰(h0(v0) − h₋(v0)),

for all c ≥ cβ⁰. Let β⁰ be such that ˆw0(0) = β⁰(h0(v0) − h₋(v0)). Then for all c > cβ⁰

ˆ

w0(c) < ˆw0(0). (4.27)

w

u

x∗0 (h0(v0), 0)

L

S0 Tβ0

λ− 0

Figure 4.7: Situation sketch for region S⁰.

Combining (4.26) and (4.27) we deduce that there exists a c > 0 such that ˆ

w0(c) < ˆw1(c). (4.28)

Lemma 6. ˆw0: [0, ∞) → R is strictly decreasing and ˆw1: [0, ∞) → R is strictly increasing.

Proof. Consider equation (4.12), subject to the conditions lim

ξ→−∞(u, w) = (h₊(v₀), 0), lim

ξ→∞(u, w) = (h₋(v₀), 0).

The eigenvalues of the linearisation around the steady states (h±(v₀), 0) are given in (4.14) and (4.15).

We will only prove that ˆw1 is strictly increasing. With a similar proof we can show that c 7→ ˆw0(c)is monotonically decreasing. For this we can rescale ξ by putting τ = −ξ, so that system (4.12) transforms into

du dτ = −w,

dw

dτ = cw + f (u, v0), (4.29)

with

lim (u, w) = (h_±(v0), 0).

(25)

In this case, the eigenvalues for (h±(v₀), 0)are

ˆλ^±₀ = c ±pc²− 4f⁰(h−(v0), v0)

2 , (4.30)

ˆλ^±₁ = c ±pc²− 4f⁰(h₊(v₀), v₀)

2 . (4.31)

Here, the stable manifold corresponds to ˆλ⁻1 and the unstable manifold corresponds to ˆλ⁺0. c 7→ ˆw0(c) will be monotonically increasing under the scaling τ = −ξ and thus monotonically decreasing for ξ.

Let c2> c1> 0and u0close to h+(v0)xed.

Remark 1. We will write Wc^u_i for the unstable manifold W^u corresponding to the dynamical system (4.12) with velocity c = ci, i = 1, 2.

Moreover we will write wi(u₀)for the solution w corresponding to a value u0

close to h+(v₀)and speed ci.

If the initial condition (u0, w0)is on Wc^u_i(h+(v0), 0), then the corresponding solution to (4.12) satises

lim

ξ→−∞

w⁰(ξ)

u⁰(ξ) = λ⁺₁(ci).

The derivative of λ⁺1(c)is dλ⁺₁

dc = −1 2+1

2

c

pc²− 4f⁰(h₊(v₀), v₀).

Since there is no c such that ^dλ_dc⁺¹ = 0and ^dλ_dc⁺¹ < 0for all c, λ⁺1(c)is monotonically decreasing for c > 0, see (4.25).

Therefore, for c2> c₁> 0, we have the inequality

λ⁺₁(c2) = −c2+pc²₂− 4f⁰(h+(v0), v0)) 2

< −c1+pc²₁− 4f⁰(h−(v0), v0) 2

= λ⁺₁(c1). (4.32)

Since w1(u0) = λ⁺₁(c1)(u0− h+(v0))and w2(u0) = λ⁺₁(c2)(u0− h+(v0))for u0

close to h+(v0), the following inequality holds

w2(u0) < w1(u0) < 0. (4.33) We now show that for every u ∈ (h−(v0), h+(v0)): w2(u) < w1(u) < 0. The unstable manifold of (h+(v0), 0) is the solution to (4.5). The slope of the line tangent to the manifold at (u, wi)is

dwi

du = −ci−f (u, v0) wi

.

Assume that Wcû1(x^∗₁) and Wcû2(x^∗₁) intersect in {(u, w)|h0(v0) + ε ≤ u ≤ h (v ), w < 0}. Since Wû(x^∗)is a C¹ submanifold, for i = 1, 2, and because

Waves in excitable media

T.M.J. van Zalen