Murray J.D. Mathematical Biology: I. An Introduction

Подождите немного. Документ загружается.

416 11. Reaction Diffusion, Chemotaxis, Nonlocal Mechanisms

Section 11.4. If we now include a reaction or dynamics term f (n) in (11.59) we get the

generalised reaction diffusion equation equivalent to (11.14). With the one space dimen-

sional scalar version of (11.59) and a logistic growth form for f (n), Cohen and Murray

(1981) have shown that the equation can exhibit steady state spatially inhomogeneous

solutions. Lara-Ochoa (1984) analysed their model in a two-dimensional setting and

showed that it reﬂects certain morphogenetic aspects of multicellular systems formed

by motile cells.

Exercises

1 Let p(x, t) be the probability that an organism initially at x = 0isatx after a

time t. In a random walk there is a slight bias to the right; that is, the probabilities of

moving to the right and left, α and β, are such that α −β = ε>0, where 0 <ε 1.

Show that the diffusion equation for the concentration c(x, t) = Qp(x, t),whereQ

particles are released at the origin at t = 0, is

∂c

∂t

+ V

∂c

∂x

= D

∂

∂x

where V and D are constants which you should deﬁne.

2 In a one-dimensional domain suppose insects are attracted to the origin x = 0andare

convected there by a constant velocity V . If the population pressure is approximated

by a density-dependent diffusion coefﬁcient D(n) = D

(n/n

)

,wheren is the

population density and D

, n

and m are positive constants, show that the model

equation for dispersal, in the absence of any population growth, is

∂n

∂t

=−

∂ J

∂x

∂

∂x

[V sgn(x)n]+D

∂

∂x





∂n

∂x



Show that if n → 0, ∂n/∂x → 0as|x |→∞a spatially inhomogeneous steady

state population density exists and can be represented by

n(x) = n



1 −

mV|x |



1/m

, if |x |≤

= 0, if |x | >

3 The larvae of the parasitic worm (Trichostrongylus retortaeformis) hatch from eggs

in sheep and rabbit excreta. The larvae disperse randomly on the grass and are conse-

quently eaten by sheep and rabbits. In the intestines the cycle starts again. Consider

the one-dimensional problem in which the larvae disperse with constant diffusion

and have a mortality proportional to the population. Show that n satisﬁes

∂n

∂t

= D

∂

∂x

−µn, D > 0,µ>0,

Exercises 417

where n is the larvae population. Find the population distribution at any x and t

arising from N

larvae being released at x = 0att = 0. Show that as t →∞the

population dies out.

If the larvae lay eggs at a rate proportional to the population of the larvae,

that is,

∂ E

∂t

= λn,λ>0,

where E(x, t) is the egg population density, show that in the limit as t →∞a

nonzero spatial distribution of eggs persists. [The result for E(x, t) is an integral

from which the asymptotic approximation can be found using Laplace’s method

(see, for example, Murray’s 1984 Asymptotic Analysis): the result gives E(x, t) ∼

O(exp[−(µ/D)

1/2

|x |]) as t →∞.]

4 Consider the density-dependent diffusion model for insect dispersal which includes

a linear death process which results in the following equation for the population

n(x, t),

∂n

∂t

= D

∂

∂x





∂n

∂x



−µn, D

> 0,µ>0.

If Q insects are released at x = 0att = 0, that is, n(x, 0) = Qδ(x), show, using

appropriate transformations in n and t, that the equation can be reduced to an equiva-

lent equation with µ = 0. Hence show that the population wavefront reaches a ﬁnite

distance x

max

from x = 0ast →∞,where

max

(µmτ

)

1/(m+2)

where

Q(

)

1/2

(

+1)

,τ

(m +2)

where  is the Gamma function.

12. Oscillator-Generated Wave Phenomena and

Central Pattern Generators

In Chapter 11 we saw how diffusion, chemotaxis and convection mechanisms could

generate spatial patterns; in Volume II we discuss mechanisms of biological pattern for-

mation extensively. In Chapter 13, and Chapter 1 and Chapter 13, Volume II we show

how diffusion effects, for example, can also generate travelling waves, which have been

used to model the spread of pest outbreaks, travelling waves of chemical concentration,

colonization of space by a population, spatial spread of epidemics and so on. The exis-

tence of such travelling waves is usually a consequence of the coupling of various effects

such as diffusion or chemotaxis or convection. There are, however, other wave phenom-

ena of a quite different kind, called kinematic waves, which exhibit wavelike spatial

patterns, which depend on the coupling of biological oscillators whose properties relat-

ing to phase or period vary spatially. The two phenomena described in this chapter are

striking, and the models we discuss are based on the experiments or biological phenom-

ena which so dramatically exhibit them. The ﬁrst involves the Belousov–Zhabotinskii

reaction and the second, which is speciﬁcally associated with the swimming of, for ex-

ample, lamprey and dogﬁsh, illustrates the very important concept of a central pattern

generator. The results we derive here apply to spatially distributed oscillators in general.

12.1 Kinematic Waves in the Belousov–Zhabotinskii Reaction

When the reactants in the oscillating Belousov reaction involve an iron catalyst (with

going to Fe

and vice versa) the oscillations are dramatically illustrated with an

appropriate dye which reﬂects the state of the catalyst: the colour change is from red

(or rather a reddish orange) to blue. When the reactants are left unstirred in a vertical

cylindrical tube horizontal bands of blue and red form. These bands usually start to

appear at the bottom of the cylinder and move slowly upwards with successive bands

moving progressively more slowly. Eventually the cylinder is ﬁlled by these bands but

with a nonuniform density, the closer to the bottom the denser the wave packing. Dif-

fusion plays a negligible role in the formation and propagation of these bands, unlike

the waves we discuss later. Beck and V

aradi (1972) provided a kinematical explanation

for these spatial patterns of bands. The analysis explaining them, which we give here,

is that of Kopell and Howard (1973). Although the analysis was originally given for the

bands observed in the Belousov reaction, and the experimental results shown in Fig-

ure 12.1(b) are also for this reaction, the phenomenon and analysis apply equally to any

12.1 Belousov–Zhabotinskii Reaction Kinematic Waves 419

biological oscillator under similar circumstances. The important point to note is that

spatial patterns can be obtained without diffusion, convection or chemotaxis playing

any role.

Consider each position in the vertical cylinder to be an independent oscillator with

period T , which may be a function of position. If these independent oscillators are

out of phase or have different frequencies then spatial patterns will appear simply as a

consequence of the spatial variation in the phase or frequency. (A simple but illustrative

physical demonstration of the phenomenon is given by a row of simple pendula all

hanging from the same horizontal rod but with a very slight gradient in their lengths.

The slight gradient in their lengths gives a slight gradient in their periods. If they are all

set swinging at the same time, then after a very short time it looks as if there is a wave

propagating along the line of pendula, the wavelength of which gets smaller and smaller

with time.)

Returning to the Belousov oscillator, the cause of a gradient in phase or frequency

can be due to a concentration gradient in one of the chemicals, or a temperature vari-

ation. The experiments (see Figure 12.1(b)) carried out by Kopell and Howard (1973)

used the former while Thoenes (1973) used the latter. The vertical chemical concentra-

tion gradient was in sulphuric acid. This resulted in a monotonic gradient in the period

of oscillation and horizontal bands appeared quite quickly, moved slowly upwards and

after a few minutes ﬁlled the cylinder. It is clear that if a barrier, impermeable to any

of the chemicals, were put in the cylinder it would neither affect the pseudowave prop-

agation nor the density of bands: spatial transport processes are simply not involved in

the generation of this spatial pattern. These ‘waves’ are indeed only pseudowaves since

nothing is actually being transported.

Let z be the spatial coordinate measured vertically from the bottom of the cylinder,

takentobez = 0, and the cylinder height to be normalised so that the top is z = 1.

Because of the initial concentration gradient there is a gradient in the oscillator period.

At position z let the period of oscillation be T (z) deﬁned for all 0 ≤ z ≤ 1. We char-

acterise the state of the oscillator by a 2π-periodic function of its phase denoted by

φ(z, t). In the Belousov reaction, for example, the front of the wave, deﬁned as the

point where φ has a speciﬁc value, can be distinguished by the sharp blue front. Let the

initial distribution of the phases be φ

(z). We can then represent the phase φ(z, t) by

φ(z, t) = ψ(z, t) +φ

(z), ψ(z, 0) = 0, (12.1)

where ψ(z, t) is a function which increases by 2π if the time t increases by the periodic

time T (z);thatis,

φ(z, t + T (z)) = 2π +ψ(z, t) +φ

(z).

Let us now take some reference phase point say, φ = 0. Deﬁne t

∗

(z) as the time at

position z at which the phase is zero; that is, it satisﬁes

0 = φ(z, t

∗

(z)) = ψ(z, t

∗

(z)) +φ

(z). (12.2)

Then, for any integer n and time t = t

∗

(z) + nT(z) we have

420 12. Oscillator-Generated Wave Phenomena

φ(z, t

∗

+nT) = ψ(z, t

∗

+nT) + φ

(z)

= 2nπ +ψ(z, t

∗

) +φ

(z)

= 2nπ,

(12.3)

using (12.1) and (12.2). So, in the (z, t) plane the point (z, t) which corresponds to the

phase 2nπ moves on the curve given by

t = t

∗

(z) + nT(z). (12.4)

We can continue the analysis with complete generality but it is just as instructive

and easier to see what is going on if we are more speciﬁc and choose T (z) tobe,say,a

smooth monotonic increasing function of z in 0 ≤ z ≤ 1. Also for simplicity, let us take

the initial distribution of phases to be a constant, which we can take to be zero; that is,

(z) = 0. This means that at t = 0 all the oscillators are in phase. From the deﬁnition

of t

∗

(z) in (12.2) this means that t

∗

(z) = 0.

If we now deﬁne

(z) = nT(z) (12.5)

then (12.3) gives

φ(z, t

(z)) = 2nπ. (12.6)

This means that t

(z) is the time at which the nth wavefront passes the point z in the

cylinder. The velocity v

(z) of this nth wavefront is given by the rate of change of the

position of the front. That is, using the last two equations,

(z) =





φ=2nπ



(z)



−1



(z)

. (12.7)

With T (z) a monotonic increasing function of z, T



(z)>0andsothenth wave-

front, the leading edge say, starts at z = 0 at time t = nT(0) and, from (12.7), propa-

gates up the cylinder at 1/n times the velocity of the ﬁrst wave. This nth wave reaches

z = 1, the top of the cylinder, at time t = nT(1). Since the T (z) we have taken here is

a monotonic increasing function of z there will be more and more waves in 0 ≤ z ≤ 1

as time goes on, since, with the velocity decreasing proportionally to 1/n from (12.7),

more waves enter at z = 0thanleavez = 1 in the same time interval.

Let us consider a speciﬁc example and take φ

(z) = 0andT (z) = 1 + z.Sothe

phase 2nπ moves, in the (z, t) plane, on the lines t = n(1 +z), and the phase φ(z, t) =

2πt/T (z) = 2πt/(1 + z). From (12.7) the velocity of the nth wave is v

= 1/n.The

space-time picture of the wavefronts, given by φ = 2nπ, is illustrated in Figure 12.1(a).

From the ﬁgure we see that at time t = 1thewaveφ = 2π enters the cylinder at z = 0

and moves up with a velocity v

= 1. At t = 2 the wave with phase 4π enters the

cylinder at z = 0: it moves with velocity v

= 1/2. At t = 3 the wave with φ = 6π

moves with velocity v

= 1/3 and so on. The ﬁrst wave takes a time t = 1 to traverse

the cylinder, the second takes a time t = 2 and so on. It is clear that as time goes on

12.1 Belousov–Zhabotinskii Reaction Kinematic Waves 421

Figure 12.1. (a) Wavefronts for the period distribution T (z) = 1+z and initial phase distribution φ

(z) = 0.

Note how many waves there are for different t’s: for t = 1.5 there is 1 wave in the cylinder, while for t = 7.5

there are 4 waves. (b) The experimental space–time situation equivalent to the theoretical results in (a). The

experiments were carried out for a Belousov reaction with an initial sulphuric acid gradient in the cylinder.

The total time in the ﬁgure is approximately 7 minutes and the vertical height about 20 cm; the ﬁgure is a

sketch from a negative. (After Kopell and Howard 1973)

more and more waves are in 0 ≤ z ≤ 1. This is clear from Figure 12.1(a) and the

experimental counterpart given in Figure 12.1(b): for example, at t = 3.5therearetwo

waves while at t = 7.5 there are 4. From the ﬁgure, it is also clear that as time increases

the waves are progressively more tightly packed nearer the bottom, z = 0.

Suppose the initial phases φ

(z) = 0; then, t

∗

(z) = 0 and from (12.4), T (z) =

[t − t

∗

(z)]/n. So, asymptotically for large time t, T (z) ∼ t/n, and the above analysis

for the illustrative example in which φ

(z) = 0 still applies asymptotically in time. Since

it is unlikely that the experimental arrangement which gave rise to Figure 12.1(b) had

a strictly linear T (z), such as used in the analysis for Figure 12.1(a), the experimental

results illustrate this asymptotic result quite dramatically.

The biological implications of the above analysis are of considerable importance.

Since biological oscillators and biological clocks are common, a time-varying spatial

422 12. Oscillator-Generated Wave Phenomena

pattern may be a consequence of a spatial variation in the oscillator parameters and not,

as might be supposed, a consequence of some reaction diffusion situation or other such

pattern formation mechanism such as we shall consider in detail in later chapters.

In this section the wave pattern is a continuously changing one. In the following

sections we investigate the possibility of a more coherent wave pattern generator of

considerable importance.

12.2 Central Pattern Generator: Experimental Facts in the

Swimming of Fish

A ﬁsh propels itself through water by a sequence of travelling waves which progress

down the ﬁsh’s body from head to tail and its speed is a function of the wave frequency.

It is the network of neurons arrayed down the back that controls the muscle movements

which generate the actual waves and coordinate them to produce the right effect. It is

a widely held lay belief that in mammals the generator, or rather the controlling nerve

centre for the rhythmic control of these waves, is the brain. However, in many animals

swimming occurs after the spinal cord has been severed from the brain—the technical

term is spinal transection. In the case of the dogﬁsh, for example, the phenomenon has

been known since the end of the 19th century. The swimming movement observed in

such situations shows the proper intersegmental muscle coordination.

The basis for the required rhythmic behaviour and its intersegmental coordination is

a central network of neurons in the spinal cord. It is known that there are neural networks

which can generate temporal sequences of signals, which here produce the required

cyclic patterns of muscle activity. Such networks are called central pattern generators

and by deﬁnition require no external input control for them to produce the required

rhythmic output. It is obvious how important it is to understand such neural control

of locomotion. However to do so requires modelling realism and, at the very least,

detailed information from experiments. The recent book on neural control of rhythmic

movements, edited by Cohen et al. (1988), is speciﬁcally about the subject matter of this

section and the theory and modelling chapters by Kopell (1988) and Rand et al. (1988)

are particularly relevant.

In the case of higher vertebrates there are possibly millions of neurons involved. So

it is clear that experimentation and its associated modelling should at least start with as

small a spinal cord as possible but one which still exhibits this posttransection activity.

Such neural activity, which produces essentially normal swimming, is called ‘ﬁctive

swimming’ or ‘ﬁctive locomotion.’ This description also includes the situation where,

even when the muscles which produce the actual locomotion are removed, the neural

output from the spinal cord is the same as that of an intact swimming ﬁsh.

Grillner (1974), Grillner and Kashin (1976) and Grillner and Wall

en (1982) present

good experimental data on the dogﬁsh. Kopell (1988) uses this work as a case study for

the theory described in detail in her paper. The lamprey, which is rather a primitive

vertebrate, was the animal used in a series of interesting and illuminating experiments

by Cohen and Wall

en (1980) (see also Cohen and Harris-Warrick 1984 and references

given there). They studied a speciﬁc species of the lamprey which varies from about

13–30 cm in length and has a spinal cord about 0.3 mm thick and 1.5 mm wide. Its

12.2 Central Pattern Generator 423

Figure 12.2. (a) Schematic diagram of the exposed lamprey spinal cord and the experimental arrangement.

The preparation was pinned through membrane tissue as shown in the schematic cross-section in (b). There

are ventral roots (VR) at each segment and the electrical activity of these was measured by the electrodes.

(From Cohen and Wall

en 1980)

advantages (or disadvantages from the lamprey’s point of view) are that it has relatively

few cells, but still has the necessary basic vertebrate organisation, and it exhibits ﬁctive

swimming behaviour. In the experiments of Cohen and Wall

en (1980) and Cohen and

Harris-Warrick (1984), dissected lengths of spinal cord from about 25–50 segments

were used: the lamprey has about 100. The animal was decapitated and the spinal cord

exposed but with most of the musculature intact. Motoneuron activity was monitored

with electrodes placed on two opposing ventral roots of a single segment; Figure 12.2

shows schematically the experimental setup of the spinal cord.

In these experiments the cord was placed in a saline solution and the ﬁctive swim-

ming, that is, the periodic rhythmic activity, was induced chemically (by L-DOPA or by

the amino acid D-glutamate); the ﬁctive swimming can go on for hours. The ventral root

(VR) recordings from the electrodes showed alternating bursts of impulses between the

left and right VR of a single segment. That is, the periodic bursts on either side of a

segment are 180

◦

out of phase. Figure 12.3 illustrates the VR activity obtained from the

left and right sides of the ventral roots from two different segment levels in the cord. An

Figure 12.3. Bursting activity recorded from the left (L) and right (R) sides of the ventral roots (VR) at

segments 7 and 19 as measured from the head-end of the specimen, which consisted of 27 segments. (From

Cohen and Wall

en 1980) The time between bursts is approximately 1 sec. Note the approximately constant

phase lag as you go from segment 7 to 19.

424 12. Oscillator-Generated Wave Phenomena

Figure 12.4. Typical swimming pattern illustrating

a propagating wave generated by a ventral root

output such as illustrated in Figure 12.3.

important point to note, and which we use in the model, is that the left and right VR of

a segment are like individual oscillators which are phase locked 180

◦

out of phase. This

intrasegmental coupling is very stable.

The results in Figure 12.3 are from an isolated piece of spinal cord consisting of

27 segments and with the numbering starting from the head side; the recordings were

taken at segments 7 and 19. The period of the bursts of activity was about 1 per second.

Another point to note is the nearly constant phase lag between the two segments: the lag

between the right VR of segment 7 and segment 19 is to a ﬁrst approximation the same

as between the left VR of segments 7 and 19. In Section 12.4 we incorporate this type of

behaviour in one of the speciﬁc cases in the model developed below. A piece of spinal

cord of only about 10 segments can produce a stable neural ﬁctive swimming output.

This periodic activity of the isolated lamprey spinal cord, which is directly related

to the undulatory wavelike movements of the swimming ﬁsh, schematically shown in

Figure 12.4, is the phenomenon controlled by the central pattern generator, which we

now wish to model. Various models can be suggested for the generation of these pat-

terns; see, for example, Kopell (1988). In the following Sections 12.3 and 12.4 we de-

scribe in detail a model proposed by Cohen et al. (1982); it has been used with consid-

erable success to explain certain experimental results associated with selective surgical

lesions in the spinal cord (Cohen and Harris-Warrick 1984).

12.3 Mathematical Model for the Central Pattern Generator

The basic characteristics of the phenomenon, as exempliﬁed by the experiments and

indicated in particular by Figure 12.3, are that the left and right VR of a segment are

phase locked oscillators and that there is approximately a constant phase lag from the

head to tail of the spinal cord. The key assumptions in the model are: (i) Each segment

in the back has associated with it a pair of neuronal oscillators each of which exhibits,

in isolation, a stable limit cycle periodic oscillation. The amplitude of the oscillation

depends only on internal parameters, and is not usually affected by external factors

such as drugs or electrical stimulations. (ii) Each of the oscillators is coupled to its

nearest neighbour but with the possibility of long range coupling; there is experimental

evidence for the latter in Buchanan and Cohen (1982).

12.3 Mathematical Model for the Central Pattern Generator 425

We saw in Chapter 7 on biological oscillators that many biochemical reaction sys-

tems can exhibit stable limit cycle periodic oscillations. It is such a biological oscillator,

or one coupled to some neuronal electrical property, which we envisage to be the driving

force in each of the oscillators associated with the spinal segments. It is not necessary

to know the actual details of the biological oscillator for our model here—we do not

know what it is in fact. As a preliminary to studying the intersegmental linking of the

oscillators we ﬁrst consider a single oscillator to set up the mathematical treatment and

notation and introduce the analytical procedure we use.

Single Oscillator and Oscillator Pair

Denote the vector of limit cycle variables of relevance by the vector

x(t) = (x

(t), x

(t),... ,x

(t)). (12.8)

Again we do not know what quantity it is that oscillates, only that something does

which gives rise to the periodic VR neuronal activity which is observed experimentally.

For example, x(t) could include the level of the neurotransmitter substance and the peri-

odically varying electric potential. We denote the vector differential equation governing

the limit cycles by

= f(x), (12.9)

where t denotes time. Just as we do not need to know the speciﬁc biological oscillator

involved, we do not require the detailed functional form of the function f(x).Letus

consider the limit cycle to be the closed orbit γ in the phase space. Using this curve

as the local coordinate system we can think of the periodic limit cycle as having a

phase θ which goes from 0 to 2π as we make a complete circuit round the closed

orbit γ . Assume that at some point P on the closed curve the bursting, which is observed

experimentally by the electrodes, occurs at the phase θ = 0. Starting at P, the phase

increases from 0 and reaches 2π when we get back to P, where bursting again occurs.

Let us further assume that the coordinate system for the limit cycle is chosen so that the

speed of the solution round γ , as measured now by dθ/dt, is constant.

The above idea of representing a limit cycle in terms of the phase can be illus-

trated by the following, admittedly contrived, simple pedagogical example. Consider

the differential equation system given by

= x

(1 − ρ) − ωx

= x

(1 − ρ) +ωx

ρ = (x

+ x

)

1/2

(12.10)

where ω is a positive constant. Although the solution of this system can be obtained triv-

ially, as we see below, a formal phase plane analysis (see Appendix A) of (12.10) shows

that (0, 0) is the only singular point and it is an unstable spiral, spiralling anticlockwise.

A conﬁned set can be found (just take ρ large), so by the Poincar

e–Bendixson theorem

a limit cycle periodic solution exists and is represented by a closed orbit in the (x

, x

)