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

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

256 7. Biological Oscillators and Switches

and sketch the domain. Indicate how the domain for periodic solutions changes as

the inhibition parameter K varies.

7 The two-variable FitzHugh–Nagumo model for space-clamped nerve axon ﬁring

with an external applied current I

= v(a −v)(v − 1) −w + I

= bv −γw,

where 0 < a < 1andb, γ and I

are positive constants. Here v is directly related

to the transmembrane potential and w is the variable which represents the effects of

the various chemical ion-generated potentials.

Determine the local maximum and minimum for the v null cline in terms

of a and I

and hence give the corresponding piecewise linear approximate form.

Show that there is a conﬁned set for the model system. Using the piece-

wise linear model, determine the conditions on the parameters such that the positive

steady state is stable but excitable. Find the conditions on the parameters, and the rel-

evant window (I

, I

) of applied currents, for the positive steady state to be linearly

unstable and hence for a limit cycle solution to exist. For a ﬁxed set of parameters

a, b and γ , ﬁnd the period of the small amplitude limit cycle when I

is just greater

than the bifurcation value I

. [Use the Hopf bifurcation result near bifurcation.]

8 Consider a simpliﬁed model for the control of testosterone secretion given by

= f (T ) − b

= b

R(t −τ) − b

where R denotes the luteinising hormone releasing hormone (LHRH), T denotes

the hormone testosterone and f (T ) is a positive monotonic decreasing function of

T . The delay τ is associated with the blood circulation time in the body and b

, b

and b

are positive constants. When τ = 0, show that the steady state is stable. Using

the method in Section 7.6 investigate the possibility of periodic solution behaviour

when τ>0.

8. BZ Oscillating Reactions

8.1 Belousov Reaction and the Field–K

os–Noyes (FKN) Model

The reaction known as the Belousov–Zhabotinskii reaction is an important oscillating

reaction discovered by the Russian Boris Belousov (1951), a biochemist, and is de-

scribed in an unpublished paper, which was contemptuously rejected by a journal editor;

at the time the accepted dogma was that oscillating reactions were simply not possible.

A translation of the original article is given in the book edited by Field and Burger

(1985). Eventually Belousov (1959) published a brief note in the obscure proceedings

of a Russian medical meeting. Basically he found oscillations in the ratio of concen-

trations of the catalyst; in Belousov’s reaction it was cerium in the oxidation of citric

acid by bromate. The oscillation manifested itself via a colour change as the cerium

changed from Ce

to Ce

although it is more dramatic with an iron ion, ferroin

where the colour is brick red when in the Fe

state and bright blue in the Fe

state.

The study of this reaction was continued by Zhabotinskii (1964) and is now known

as the Belousov–Zhabotinskii reaction or simply the BZ reaction. When the details of

this important reaction and some of its dramatic oscillatory and wavelike properties

ﬁnally reached the West in the 1970’s it provoked widespread interest and research.

Belousov’s seminal work was ﬁnally, but posthumously, recognised in 1980 by his be-

ing awarded the Lenin Prize. Winfree (1984) gives a brief interesting description of

the history of the Belousov–Zhabotinskii reaction. When the reactants can also diffuse

a diverse menagerie of complex patterns can be formed and it is the latter which has

sustained the continuing widespread interest among both biological and physical scien-

tists. It provided an enormous incentive to those who were interested in the fundamental

question of spatiotemporal self-organisation with all the implications for the generation

of biological pattern and form. We discuss such aspects in considerable detail in later

chapters.

Before discussing the reaction in detail it is highly pertinent to mention its analogy

with real biological oscillators, not at the molecular level, of course, but in the similarity

it has with an increasing number of real biological cyclical behaviour. Approximate

equations for the BZ reaction are identical to those that have arisen in a realistic model

for the periodic signalling exhibited by cells during the self-organisation of the slime

mould Dictyostelium discoideum (see, for example, the book by Goldbeter 1996). There

is, of course, no analogy at the molecular level but nevertheless a knowledge of the

behaviour from similar equations can still be very useful. Winfree (1987) also used it in a

258 8. BZ Oscillating Reactions

highly original and elegant way in modelling three-dimensional activity in the ventricle

of the heart. The analogy has also been applied by Tyson (1991) in his investigations

into the molecular biology of the cell cycle. This reaction is certainly not just some

academic curiosity.

There are now many such chemical reactions which can exhibit periodic behaviour

and the term BZ reaction now refers to a general class of such reactions, essentially

where an organic substance is oxidised by bromate ions facilitated by a metal ion in an

acid medium. Typical metal ions are cerium and ferroin. Although it is a chemical rather

than a biological oscillator the BZ reaction is now considered the prototype oscillator.

The detailed reactions involved are more or less understood, as are many, but certainly

not all, of the complex spatial phenomena it can exhibit: we shall describe some of the

wavelike properties in Chapter 1, Volume II. The book of articles edited by Field and

Burger (1985) is a good and varied introduction to the BZ reaction. It also has several

articles on chemical oscillators and wave phenomena. The book by Winfree (2000),

among other things, also discusses some of the reaction’s properties, both temporal

and spatial. The article by Tyson (1994) is a succinct and more recent summary of the

phenomena exhibited by the BZ kinetics, both temporal and spatial. In this chapter we

consider the reaction in some detail not only because of its seminal importance in the

ﬁeld, but also because it illustrates techniques of analysis which have wide applicability.

References to the more detailed kinetics are given at the appropriate places. Almost all

the phenomena theoretically exhibited by reaction and reaction diffusion mechanisms

have been found in this real and practical reaction—but many of these only after the

mathematics predicted them.

The BZ reaction (the general class of them) is probably the most widely studied

oscillating reaction both theoretically and experimentally. Here we brieﬂy describe the

key steps in the reaction and develop the Field–K

os–Noyes (Field and Noyes 1974)

model system which quantitatively mimics the actual chemical reactions (Field et al.

1972). The models for the BZ reaction are prototypes to study since the theoretical

developments can be tested against experiments. The experience gained from this is

directly transferable to biochemical oscillators as mentioned above. The literature on

the subject is now large. A succinct review of the detailed reaction and its properties is

given by Tyson (1994).

In the original Belousov (1951) reaction, the basic mechanism consists of the ox-

idation of malonic acid, in an acid medium, by bromate ions, BrO

−

, and catalyzed

by cerium, which has two states Ce

and Ce

. Sustained periodic oscillations are

observed in the cerium ions. With other metal ion catalysts and appropriate dyes, for

example, iron Fe

and Fe

and phenanthroline, the regular periodic colour change

is visually dramatic, oscillating as mentioned between red and blue. It is not only the

catalyst ion concentrations which vary with time, of course, other reactants also vary.

Figure 8.1 illustrates the temporal variations in the bromide ion concentration [Br

−

]and

the cerium ion concentration ratio [Ce

]/[Ce

] measured by Field et al. (1972), who

studied the mechanism in depth; see Tyson (1994) for more references and technical

details.

Basically the reaction can be separated into two parts, say I and II, and the concen-

tration [Br

−

] determines which is dominant at any time. When [Br

−

] is high, near A in

Figure 8.1, I is dominant and during this stage Br

−

is consumed; that is, we move along

8.1 Belousov Reaction and Field–K

os–Noyes Model 259

Figure 8.1. Experimentally measured periodic limit cycle type of temporal variation in the concentrations in

the ratio of the cerium metal ion concentration [Ce

]/[Ce

] and the bromide ion concentration [Br

−

]in

the Belousov–Zhabotinskii reaction. (Redrawn from Field et al. 1972)

AB, and the cerium ion is mainly in the Ce

state. As [Br

−

] decreases further it passes

through a critical value, B, and then drops quickly to a low level, that is, C in Figure 8.1.

At this stage process II takes over from I. During II the Ce

changes to Ce

.However,

in the II-process Ce

reacts to produce Br

−

again while it reverts to the Ce

state.

Now [Br

−

] increases, that is, along CDA, and, when its value is sufﬁciently high, pro-

cess I again becomes dominant. The whole sequence is continually repeated and hence

produces the observed oscillations. This qualitative description is not, of course, sufﬁ-

cient to show that the reaction will actually oscillate. The two pathways, I and II, could

simply reach some steady state of coexistence. In fact for certain parameter ranges this

is exactly what can happen. The parameter range required for oscillatory behaviour is

derived in detail below and depends on a careful analysis of the kinetics equations.

The rapid variation along BC and DAis typical of a relaxation oscillator, which is

just an oscillator in which parts of the limit cycle are traversed quickly. This behaviour

suggests a particular asymptotic technique which often allows us to get analytical results

for the period in terms of the parameters; we discuss this later in Section 8.4.

Although there are many reactions involved they can be rationally reduced to 5 key

reactions, with known values for the rate constants, which capture the basic elements of

the mechanism. These ﬁve reactions can then be represented by a 3-chemical system in

which the overall rate constants can be assigned with reasonable conﬁdence. This model

is known as the Field–K

os–Noyes or FKN model and is the speciﬁc model proposed

by Field and Noyes (1974) based on the Field–K

os–Noyes (1972) mechanism. We

give this simpler model system and derive the 3-species model. A complete derivation

from the chemistry together with estimates for the various rate constants are given by

Tyson (1994).

The key chemical elements in the 5-reaction FKN model are

X = HBrO

, Y = Br

−

, Z = Ce

A = BrO

−

, P = HOBr,

(8.1)

and the model reactions can be approximated by the sequence

260 8. BZ Oscillating Reactions

A +Y

→ X + P, X + Y

→ 2P

A + X

→ 2X + 2Z, 2X

→ A + PZ

→ fY,

(8.2)

where the rate constants k

,... ,k

are known and f is a stoichiometric factor, usually

taken to be 0.5. The ﬁrst two reactions are roughly equivalent to the process I, described

above, while the last three relate approximately to process II. The form used by Field

and Noyes (1974) is actually slightly different in that the last reaction is B + Z →

fY + Q where B represents organic acids and Q is CO

. [B] is a constant and can

be incorporated into the rate constant in the model equation system. It is reasonable to

take the concentration [A] of the bromate ion to be constant; the concentration [P] is

not of interest here. So, using the Law of Mass Action on (8.2), we get the following

third-order system of kinetics equations for the concentrations, denoted by lowercase

letters:

= k

ay − k

xy +k

ax − k

=−k

ay − k

xy + fk

= 2k

ax − k

(8.3)

This oscillator system is sometimes referred to as the ‘Oregonator’ since it exhibits limit

cycle oscillations and the research by Field et al. (1972) was done at the University of

Oregon.

Oscillatory behaviour of (8.3) depends critically on the parameters involved. For

example, if k

= 0, the bromide ion (Br

−

) concentration y decays to zero according to

the second equation, so no oscillations can occur. On the other hand if f = 0.5andk

very large, the last reaction in (8.2) is very fast and the third and ﬁfth reactions in (8.2)

effectively collapse into the single reaction

A + X

→ 2X +Y.

The system then reduces to a two-species mechanism, which is bimolecular, and so

cannot oscillate (Hanusse 1972). There is clearly a domain in the ( f, k

) plane where

periodic behaviour is not possible.

The only sensible way to analyse the system (8.3) is in a dimensionless form. There

are usually several ways to nondimensionalise the equations; see, for example, Murray

(1977) and Tyson (1982, 1985), both of whom give a fuller description of the chemistry

and justiﬁcation for the model. Here we give only the more recent one suggested by

Tyson (1985) which incorporates the estimated values of the rate constants in the above

model (8.3). Often different nondimensionalisations highlight different features of the

oscillator. Following Tyson (1985), introduce

∗

, y

∗

, z

∗

, t

∗

8.2 FKN Model and Limit Cycle Solutions 261

≈ 1.2 ×10

−7

M, y

≈ 6 ×10

−7

2(k

≈ 5 ×10

−3

M, t

≈ 50s,

ε =

≈ 5 ×10

−5

,δ=

≈ 2 ×10

−4

q =

≈ 8 ×10

−4

,(f ≈ 0.5)

(8.4)

which we substitute into (8.3). Field and Noyes (1974) suggested the value for f , based

on experiment. As said before, since the model telescopes a number of reactions the

parameters cannot be given unequivocally; the values are ‘best’ estimates. For our pur-

poses here we need only the fact that ε, δ and q are small. With (8.4) we get the fol-

lowing dimensionless system, where for algebraic convenience we have omitted the

asterisks,

= qy − xy + x(1 − x),

=−qy − xy + 2 fz,

= x − z.

(8.5)

In vector form with r = (x, y, z)

we can write this as

= F(r;ε, δ, q, f ) =







−1

(qy − xy + x − x

)

−1

(−qy − xy +2 fz)

x − z







. (8.6)

8.2 Linear Stability Analysis of the FKN Model and Existence of

Limit Cycle Solutions

Even though the system (8.5) is third-order, the linear stability analysis procedure is

standard and described in detail in Chapter 3; namely, ﬁrst ﬁnd the positive steady state

or states, determine the eigenvalues of the linear stability matrix and look for a conﬁned

set, which is a ﬁnite closed surface S enclosing the steady state such that any solution at

time t

which lies inside S always remains there for all t > t

. This was done by Murray

(1974) for the original, but only slightly different, FKN equation system; it is his type

of analysis we apply here to (8.6).

The nonnegative steady states (x

, y

, z

) of (8.5) are given by setting the left-hand

sides to zero and solving the resulting system of algebraic equations, to get

262 8. BZ Oscillating Reactions

(0, 0, 0) or z

= x

, y

2 fx

q + x



(1 − 2 f − q) +[(1 −2 f −q)

+4q(1 + 2 f )]

1/2



(8.7)

The other nonzero steady state is negative.

Linearising about (0, 0, 0) we obtain the stability matrix A with eigenvalues λ given

| A −λI |=



−1

−λ qε

−1

0 −qδ

−1

−λ 2 f δ

−1

10−1 −λ



= 0

⇒ λ

+λ

(1 +qδ

−1

−ε

−1

) − λ[ε

−1

(1 + qδ

−1

) − qδ

−1

]−

q(1 + 2 f )

εδ

= 0.

If we simply sketch the left-hand side of this cubic as a function of λ for λ ≥ 0wesee

that there is at least one positive root. Alternatively note that the product of the three

roots is q(1+2 f )/εδ > 0, which implies the same thing. Thus the steady state (0, 0, 0)

is always linearly unstable.

If we now linearise (8.5) about the positive steady state (x

, y

, z

) in (8.7) the

eigenvalues λ of its stability matrix are given, after a little algebra, by

| A −λI |=



1 −2x

− y

−λ

q − x

−y

−

−λ

2 f

10−1 −λ



= 0

⇒ λ

+ Aλ

+ Bλ +C = 0,

(8.8)

where, on using the quadratic for x

, the simultaneous equations for x

, y

from (8.6)

and some tedious but elementary algebra, we get

A = 1 +

q + x

E = 2x

+ y

−1 =

+q(x

+2 f )

q + x

> 0,

B =

q + x

(q + x

)E + y

(q − x

)

εδ

C =

(q + x

)E − 2 f (q − x

) + y

(q − x

)

εδ

+q(2 f + 1)

εδ

> 0.

(8.9)

Note that A > 0, since E > 0, and that C > 0, on using the expression for x

from (8.7);

8.2 FKN Model and Limit Cycle Solutions 263

B can be positive or negative. It follows from Descartes’ rule of signs (see Appendix B)

that at least one eigenvalue λ in (8.8) is real and negative. The remaining necessary and

sufﬁcient condition for all of the solutions λ to have negative real parts is, from the

Routh–Hurwitz conditions (see Appendix B), AB − C > 0. Substituting for A, B and

C from (8.9) gives a quadratic in 1/δ for the left-hand side and hence the condition for

stability of the positive steady state in (8.7). This is given by

AB −C = φ(δ, f,ε) =

Nδ

+ Mδ + L

> 0,

L = (q + x

)



(q + x

) +

(1 −q − 4 f ) + 2q(1 +3 f )



N =[x

+q(x

+2 f )]

1 + E/ε

ε(q + x

)

> 0

(8.10)

with M also determined as a function of x

, f , q and ε; we do not require it in the

subsequent analysis and so do not give it here. With x

from (8.7), L, M and N are

functions of f , q and ε. Thus for the steady state to be linearly unstable, δ, f and ε must

lie in a domain in (δ, f,ε) space where φ(δ, f,ε) < 0. The boundary or bifurcation

surface in (δ, f,ε)space is given by φ(δ, f,ε) = 0.

We can get an indication of the eigenvalue behaviour asymptotically for large pos-

itive and negative B.IfB  1 the asymptotic solutions of (8.8) are given by

λ ∼−

, −

±i

√

B, (8.11)

while if B < 0and| B |1,

λ ∼

| B |

, ±



| B |. (8.12)

So, for large positive B condition (8.10) is satisﬁed and from (8.11), Re λ<0andthe

steady state is linearly stable, while if B is large and negative, it is unstable.

When the parameters are such that B = C/A, the bifurcation situation, we can

solve for the roots λ in (8.8), namely,

λ =−A, ±i

√

B, when B =

. (8.13)

If B = (C/A) − ω,0<ω 1, it can be seen by looking for asymptotic solutions to

(8.8) in the form λ =±i(C/A)

1/2

+ O(ω) that the O(ω) term has a positive real part.

Thus, near the bifurcation surface in the unstable region, the steady state is unstable by

growing oscillations. The conditions of the Hopf bifurcation theorem (see, for example,

Strogatz 1994) are satisﬁed and so in the vicinity of the surface φ(δ, f,ε) = 0the

system exhibits a small amplitude limit cycle solution with period

T =

2π





1/2

. (8.14)

264 8. BZ Oscillating Reactions

With the parameter values obtained from experiment, and given in (8.4), the amplitudes

of the oscillations in fact are not small, so the last expression is of pedagogical rather

than practical use. However, the bifurcation surface φ(δ, f,ε) = 0 given by (8.10) is of

practical use and this we must discuss further.

The surface φ(δ, f,ε) = 0, namely, Nδ

+ Mδ + L = 0 in (8.10), where N > 0,

is quadratic in δ. So, for the steady state to be unstable, that is, δ, f and ε make φ<0,

δ must be such that

0 <δ<



−M +[M

−4LN]

1/2



. (8.15)

But there is a nonzero range of positive δ only if the right-hand side of this inequality is

positive, which requires, with L from (8.10),

L = (q + x

)



(q + x

) +

(1 −q − 4 f ) + 2q(1 +3 f )



< 0.

With x

from (8.7) this gives an algebraic equation relating f , q and ε. From (8.4) ε is

small, so to ﬁrst-order the last inequality gives the numerator in the ε

−1

term equal to

zero, which reduces to

(1 −4 f − q){(1 −2 f −q) +[(1 −2 f − q)

+4q(1 +2 f )]

1/2

}

+4q(1 + 3 f )<0,

(8.16)

which (with an equals sign in place of <) deﬁnes the critical f , f

,foragivenq. In fact

there are two critical f ’s, namely,

and

and f must lie between them:

< f <

With q = 8 ×10

−4

we can determine accurate values for

and

by exploiting the

fact that 0 < q  1. The critical f

are given by (8.16) on replacing the inequality sign

with an equals sign. Suppose ﬁrst that (1 −2 f −q)>0; that is, 2 f < 1toO(1).Then

on letting q → 0,

(1 −4 f )(1 −2 f ) ≈ 0 ⇒

≈

With (1 −2 f −q)<0, that is, 2 f > 1toO(1), the limiting situation q → 0hastobe

done carefully; this gives

.ToO(q), (8.16), again with an equals sign, now becomes

(1 −4 f )



−(2 f + q −1) +(2 f + q −1)



1 +

2q(1 + 2 f )

(2 f + q −1)



+4q(1 + 3 f ) ≈ 0,

which reduces to

8.3 Nonlocal Stability of the FKN Model 265

4 f

−4 f − 1 ≈ 0 ⇒

≈

1 +

√

So, for small q, the range of f for which the positive steady state in (8.7) is linearly

unstable is

≈

< f <

≈

1 +

√

. (8.17)

Finally the stability bifurcation curve of δ against f for each ε is given by (8.15),

namely,

δ =



−M +[M

−4LN]

1/2



. (8.18)

with

< f <

, where the critical f

are obtained from (8.16), and L, M and N

are deﬁned by (8.10), with (8.9), in terms of f , q and ε.

8.3 Nonlocal Stability of the FKN Model

We showed in the last section that for each ε,ifδ and f lie in the appropriate domain, the

positive steady state is linearly unstable, and indeed by growing oscillations if (δ, f ) are

close to the bifurcation curve. Wherever (δ, f ) lies in the unstable domain, we must now

consider global stability. Even though we do not have the equivalent of the Poincar

e–

Bendixson theorem here since we are dealing with a third-order system, the existence

of a periodic solution with ﬁnite amplitudes requires the system to have a conﬁned set,

S say. That is, with n the unit outward normal to S,wemusthave

n ·

< 0, r on S, (8.19)

where dr/dt is given by (8.6).

Although the existence of a conﬁned set and a single unstable steady state (by grow-

ing oscillations) is not sufﬁcient to prove the existence of a periodic limit cycle, they

give sufﬁcient encouragement to pursue the analysis further. With three equations it is

possible to have chaotic solutions, such as we found with discrete models in Chapter 2.

Chaotic behaviour in this sense has been found in the BZ reaction (see, for example,

Scott 1991). (The now classical Lorenz (1963) model is a 3-equation system and it ex-

hibits chaos; it models a ﬂuid dynamics situation.) Hastings and Murray (1975) gave a

rigourous proof, together with a procedure which determines the general trajectory path

over a cycle, showing that the FKN model system (8.5) possesses at least one limit cycle

periodic solution. The procedure they developed has wider applications to a fairly broad

class of feedback control systems as has been demonstrated by Hastings et al. (1977).

Let us look for the simplest surface S, namely, a rectangular box deﬁned by the

faces