Chemical instabilities and relaxation oscillations

.T. theor. Biol. (1971) 32, 283-292

Chemical Instabilities and Relaxation Oscillations R. LAVENDA,?

G. NICOLS

AND

M.

HERSCHKOWITZ-KAUFMAN

FacttltP des Sciences de I’ UniversitP Libre de Bru.w/fes, Bruxelles. Belgium (Received 14 December 1970) The behaviour of a chemical system beyond a non-equilibrium unstable transition is analysed in the limit of a large separation of time scales characterizing the individual chemical steps. It is found that the system may exhibit in this limit relaxation, or quasi-discontinuous type oscil!ations. The properties of certain types of biological regulatory phenomena and the results of recent experiments on the Zhabotinski reaction are discussed within the framework of this result.

1. Introduction In a recent paper (Lefever & Nicolis, 1971) R. Lefever and one of the present authors discussed some of the implications of dissipative instabilities (cf. Prigogine, 1969) in the temporal organization of chemical systems. The main ideas have been illustrated by considering the following model of coupled chemical reactions A -+ X, B+X -+ Y+D, 2x+y -+ 3x, X --, E, (1) wherein the inverse reaction rates are neglected and the initial and final product concentrations A, B, D, E are maintained time and space-independent. The system described by (1) is therefore necessarily infinitely far from thermodynamic equilibrium (Prigogine, 1967). The chemical kinetic equations giving the evolution of the intermediate products read (we assume that the system is maintained spatially uniform): 8 = A-(B+ l)XfX2Y, y = BX-X”Y.

(2a) (2b)

t Present address: Department of Physical Chemistry, The Hebrew University. Jerusalem, Israel. T.R. 28 1 IV

284

B.

LAVENDA,

G.

NICOLIS

AND

M.

HERSCHKOWITZ-KAUFMAN

The analytic and numerical study of these equations has revealed the following features (Lefever & Nicolis, 1971). (i) The system admits a single steady state X, = A,

Y, = B/A.

(3)

(ii) The steady state (3) is unstable with respect to space-independent infinitesimal perturbations as long as with

B > B,, B, = A2+1.

(4)

(iii)

Beyond instability the system evolves to a new stable regime such that the concentrations of X and Y are periodic in time. These oscillations are of the limit cycle type, i.e. they represent purely non-linear effects. Their characteristics (period, amplitude) are independent of the initial conditions and only depend on the intrinsic properties of the system, i.e. the values of the initial and final product concentrations, the kinetic constants and so on. (iv) The occurrence of the chemical instabilities and the subsequent evolution of the system to a time-organized state are only possible in non-linear systems maintained in a finite distance from thermodynamic equilibrium. (v) The stability and uniqueness of periodic trajectories in some chemical systems beyond instabilities may suggest the importance of these effects in biological regulatory phenomena (cf. also Walter, 1970). The object of this note is to discuss the properties of system (1) in the limit when the second step becomes very fast. Formally this may be expressed by A = finite,

B/A --$ cn,

(5)

considering that all kinetic constants are equal to one. Relation (5) implies therefore that the time scales characterizing the individual chemical steps are largely separated.J From a purely mathematical viewpoint we expect the system to perform in this limit refuxation oscillations (Minorsky, 1962). In section 2 this is indeed shown to be the case. An asymptotic evaluation is also given of the period and amplitude of these oscillations. In section 3 we compile a few characteristic properties of the motion such as constants of motion, direction of rotation and so on. Section 4 is devoted to an interpretation of the results of the previous sections in more physical terms. It is conjectured that several 1 We notice that in thelimit (5) the instability condition (4) is always satisfied. The system therefore undergoes necessarily a stable, periodic motion.

CHEMICAL

INSTABILITIES

AND

RELAXATION

OSCILLATIONS

285

important biological phenomena may well exhibit at a certain stage a relaxation oscillation type of behaviour. Experimental evidence of the occurrence of relaxation oscillations is also discussed.

2. Relaxation Oscillations

in Model (1)

We want to study the behaviour of equations (2) in the limit (5). To this end it is convenient to introduce the new dependent variable z = *Y+Y. Equations

with

(6)

(2) and (6) then yield 2 = A+Y-Z = F(Y,Z), &jr = (Z- I’)[1 - l’(Z- Y,/B] c = B-l < 1.

= G( Y. Z),

(7a) (7b) (7c)

This system of equations may be studied by the classical methods of the theory of relaxation oscillations (Minorsky, 1962). The main point is to observe that, because of the smallness of L Y will be very large (formally Y + ~8) everywhere except on the curie G(X,%) Similarly

= 0.

the integral curves (trajectories) dZ dY

F( y, a & G(Y,Z)

(81

of the system (9)

will be practically horizontal everywhere in the (Z. Y) plane except on the curve G = 0. In Fig. I we have plotted the curve (8) for the numerical values A = 8.2, B = 77. From equation (7b) it may be seen that this curve consists of two branches Z, and Z, given by: z, = Y. zz = Y t-B/,‘.

(10)

Branch Z, has a parabolic form with Z,,,;, = 2\‘B, Y,,,i, = t.2. The steady state Y, N 9.4, Z, = 17.6 [cf. equation (3)] lies on the branch Z, and is indicated in Fig. 1 by a small circle. It is also instructive to draw the direction field on both sides of the branches. This is done by studying the sign of G in various regions of the (Z, Y) plane. The result is indicated in Fig. 1 by arrows. It is seen that near the branches Z, and Z; the direction field is directed toward these branches on both sides, while on Zs it is directed away from that branch. This shows that Z, and Zi are stable and Z; is unstable (Z; being that region of Z, where Y > Ymin).

286

B.

LAVENDA,

G.

/1 L.. _---- Ia---50 ‘;,,i,,

NICOLIS

AND

M.

100

HERSCHKOWITZ-KAUFMAN

I50

20

Y

FIG. 1. Graphical representation of curve (8) for A 7 8.2, B .= 77. Point (I’,, ZO) corresponds to the steady state solution (5). Dotted lines indicate the unstable branch of the curve.

We are now in position to follow the motion of some initial point, e.g. P, in the phase plane. As the direction field is practically horizontal, P will evolve until it meets the curve 2; at P’. It will then follow branch Z; in a prescribed direction (cf. section 3) until it reaches the point of minimum Q.t At this point the situation is as follows. Q is not a steady state position, therefore the representative point arriving at Q has to evolve further. On the other hand it cannot follow branch ZG which is unstable. Hence the only issue is to follow the (available) horizontal direction at point Q until it arrives at point R on branch Z,. The point then begins to follow branch 1 in the prescribed direction. We observe that the part QR of the path is non-analytic while the paths Z; and Z, are analytic. How long will the representative point move on Z, ? We observe that around Z,, X is practically zero. Consider now equation (2a) with X remaint Strictly speaking the representative point will not follow exact/y the branch 2; but rather a path very close to it.

CHEMICAL

INSTABILITIES

AND

RELAXATION

OSCILLATIONS

287

ing small during a certain period of time. This means that (for Y not too large) x N A-@+1)X 2: 0, i.e. x 2: A/(Bf 1). (11) Meantime Y increases and at some point the positive term X2 Y tends to increase ?i. i.e. X itself. This will happen at a point where XZY-@+1)X+.4 > 0 ( 13) for all X around the (practically zero) value (I 1). Inequality (12) implies that Y 2 @+1)2/4/l. (13) At the point on branch 2, corresponding to this value of Y the system will no longer be able to follow branch Z,. The only (non-analytic) possibility will again be a finite jump towards branch Zi. Summarizing, we obtain a piecewise analytic limit cycle QRSTQ containing the unstable steady state. It is clear that, as soon as X starts increasing Y will decrease [cf. equation (2b)]. The value of Y determined by the equality sign in equation (13) is therefore the maximum value of Y. i.e. it determines the amplitude of the oscillation of constituent Y: Yi”OX= (B+ l)li4A. (14) Let us finally estimate the period T of the motion along the limit cycle. From equation (7a) we obtain: dt = dZ/(.4 + 1’ -Z). Or, integrating over a period, dZ T=

s

TQ+QR+RS+ST

(15)

A+

y-z’

Now during the transitions QR, ST the velocity Yis practically infinite and Z remains practically constant. These transitions happen therefore during a negligible time. It follows that T = T,+T,,

__dZ ..- TI = j 4+Y-z’

i16)

TQ

T, = j ._ -!?-~m, RS .4 + Y - z

In order to calculate Tl explicitly small. This yields

we observe that on branch Z;, Y remains

288

B.

LAVENDA,

G.

NICOLIS

AND

M.

HERSCHKOWITZ-KAUFMAN

i.e.

T = ,n A-(B+1)*/4A 1 A-2&

(17a)

’

For B/A -+ co this reduces to Tl E 2 In B-

The calculation Thus

In (8A&). (17b) of T2 is based on the observation that on branch Z,, Y = Z. T*

=

zmax dZ J .A-

=

1

A

[(B + 1)2/4A - 2yiBj

LliIl

(18a)

For B/A -+ co this reduces to: T, N B2/4A2.

(18b)

The kinetic equations (2) for X and Y have also been integrated numerically for B = 77, A = 8.2. The results are shown in Figs 2 and 3. We obtain indeed a piecewise analytic limit cycle and a discontinuous oscillation with a maximum amplitude X inax = Lxx 21 189 (19a) and a period T r 253.

(19b)

On the other hand the asymptotic formulae (14), (17b) and (18b) yield: Xf& = Yzax N 185, GOa) i.e.

T,“’ N 2.3, Tas - 24.3,

T;” 21 22,

in very good agreement with the numerical results (19). Moreover one may verify from Figs 2 and 3 that the parts of motion corresponding to the quasidiscontinuous jumps [parts (b) and (d) in Fig. 21 are indeed practically instantaneous, in agreement with the arguments given earlier in this section. 3. Physical Properties

of Chemical Relaxation

Oscillations

As in all problems involving oscillations around a steady state, in the case of relaxation oscillations also the direction of rotation is determined by the chemical mechanism itself independently of the initial conditions. This is in agreement with the thermodynamic evolution criterion of Glansdorff and Prigogine (Prigogine, 1967) and is also verified by the results of the numerical computation (cf. Fig. 2 where the direction of rotation is indicated by an arrow). On inspecting Figs 1,2 and 3 one may also deduce further interesting properties of chemical relaxation oscillations. We observe that during the slow

Concn

of

y (arb~trarv

udsl

FIG. 2. Numerically computed limit cycle in the X- Y plane for A : 8.2, B -T 77. The paths (c) and (a) correspond to branches Z; and Z1 of Fig. I while the paths (d) and (b) correspond to the horizontal segments QR and ST.

r

FIG. 3. Time variation of component X as computed by numerical integration of system (2). The low concentration region corresponds to path (a) of Fig. 2 or branch Z1 of Fig. 1.

290

B.

LAVENDA,

G.

NICOLIS

AND

M.

HERSCHKOWITZ-KAUFMAh‘

parts of the motion [parts (a) and (c) in Fig. 21 the system is at a practically zero level of one chemical component and at a great excess of the other. On the other hand during the rapid motion parts [parts (b) and (d) in Fig. 21 both X and Y vary from a practically zero level to their maxima in such a way that X(t)+ Y(t) = 2 = constant. (21) Indeed parts (b) and (d) are practically straight lines with a slope of 45”. The system therefore exhibits a conservatiorl property during that part of the motion. This property has also some interesting thermodynamic implications. In a recent paper (Lefever & Nicolis, 1971) it has been shown that the average entropy production, 8 over a part of the motion is related to the steady state entropy production, co by -. .-. --. G = c-a,

=-

lnXd~~+lnY~

>

,

(22)

and that (W-r = 0, (23) where in equation (23) the average is taken over one complete period of the motion. Consider now equation (22) over the rapid parts of the motion. Using equation (21) we obtain G = d$ln(Z/Y-l)

=-$XlnX+Yln

Now along the path (d) (cf. Fig. 2) one has dY - ’ 0, dt i.e. Aa(,, < 0. Similarly, along path (b)

ii&) = d”, Y In Z-Y Y i -----ZIn(Z-Y) =

Ymin111

z-

ymin

~ymin

Y)

6’4

(25)

1

- - Y,,, In
-Z In (Z - Ymi,) + Z In (Z - Y,,,). max

Now along (b): Y,,, 2: Z. Therefore, retaining only dominant terms we obtain : LG@, = Ymin In (z/ymin) O6) > 0. From inequalities (25) and (26) we see that in the low concentration region [branch (d)] the system dissipates on the average at a lesser amount than in

CHEMICAL

INSTABILITIES

AND

RELAXATION

OSCILLATIONS

291

the steady state. In the high concentration region the situation is the opposite: the average dissipation is higher than that at the steady state. As Ymin 2: Xmir, % B- ’ jcf. equation (1 l)] we see that in this region the excess entropy production varies as Be’ In B, i.e. it attains a maximum and then goes to zero for B extremely large. In conclusion we see from the analysis of this section that at different stages ,)f the motion the system behaves quite differently as far as its functional and llissipative properties are concerned.

4. Biological Implications Relaxation oscillations are well known effects in the theory of electric circuits and in certain mechanical devices and have a great number of applicatlons in these fields (Minorsky, 1962; Andronov, Vitt & Khaikin, 1966). In this section we wish to emphasize that the physical properties of chemical relaxation oscillations described in section 3 may also be of great importance for certain types of biochemical reactions. Let us first observe that the occurrence of events with large time scale separation, which is a necessary condition for the occurrence of relaxation oscillations. is far from being exceptional in biological systems. Elementary steps involving enzyme catalysis or substrate fixation on one enzymatic form are usually very rapid compared. e.g. to nonenzymic conversion, diffusion controlled steps or conversion from one (active) macromolecular conformation to another (inactive) one. We expect therefore that a number of biochemical reactions of metabolic type could, in principle, exhibit discontinuous oscillations. On the other hand the function of a relaxation oscillator having the capacity to produce a great amount of a particular substance during a very short period of time [cf. Figs 2 and 3 and equation (21)] could well be very significant for regulatory processes involving thresholds and being sharply localized in time. On a different level, the existence of autonomous relaxation oscillators has been conjectured in a recent theory of spatial and temporal organization of developing systems (Goodwin & Cohen. 1969; Cohen, 1970). The important role of such oscillations is, again, that the oscillator could “fire” or “discharge” 11~a very sharp manner when a critical concentration is reached (Cohen. 1970). Experimentally, relaxation oscillations have been observed in an organic reaction in which cerium ions catalyse the oxidation of analogs of malonic: acid by bromate (Zaikin & Zhabotinski, 1970). When this reaction occurs in a thin layer of unstirred solution a propagation of concentration waves is observed. At every point in space the system presents relaxation oscillations

292

B.

LAVENDA,

G.

NICOLIS

AND

M.

HERSCHKOWITZ-KAUFhiAN

involving a short oxidation phase and a longer reduction phase (Zaikin Cyr Zhabotinski, 1970). The analogy with the qualitative picture given in Figs 2 and 3 is striking.

The analogy goes even further:

indeed, quite recently it has

been shown that chemical mechanisms such as (1) give rise also to sustained chemical waves in the same range of parameters for which relaxation oscillations occur (Herschkowitz-Kaufman & Nicolis, 1970). Further comments on the physical implications of relaxation oscillations and chemical waves are also made in this reference. This research has been supported, in part by the Air Force Office of Scientific Research(S.R.P.P.) through the European Office of Aerospace Research,OAR, United StatesAir Force, under grant number EOOAR 69-0058. REFERENCES A. A., VITT, A. A. PergamonPress.

ANDRONOV,

& KHAIKIN,

S.

E. (1966). Theory

qf Oscillators.

Oxford:

COHEN, M. H. (1970). “Models of Clocks and Maps in Developing Organisms”. University of Chicago preprint. GOODWIN, B. C. & COHEN, M. H. (1969). J. theor. BioL 25, 49. HERSCHKOWITZ-KAUFMAN, M. & NICOLIS, G. (1971). J. chenl. P~ZJX (in the press). LEFEVER, R. & NICOLIS, G. (1971). J. theor. Biol. 30, 267. MINORSKY, N. (1962). Non-Linear Oscillations. Princeton: van Nostrand. PRIGOGINE, I. (1967). Introduction to Thermodynanlics of Irreversible Processes, 3rd ed. New York: Interscience-Wiley. PRIG~CINE I. (1969). In Theoretical Physics and Biology. (M. Marois, ed.) Amsterdam: North Holland Publishing Co. WALTER, C. F. (1970). J. theor. Biol. 27, 259. ZAIKIN, A. N. & ZHABOTINSKI, A. M. (1970). Nature, Land. 225, 535.