Chemical oscillators: the BZ reaction

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2006

Introduction

In 1951, Belousov discovered a chemical reaction in which the concentrations of a catalyst and several reactants exhibited stable oscillations and wave phenomena, but his paper was rejected.  Although Belousov eventually published a brief report in the proceedings of an obscure Russian medical conference (1959) and Zhabotinskii continued to study this reaction in the sixties and early seventies, its importance was not widely recognized until western researchers, such as Field and Noyes, began to investigate its spatial and temporal properties in greater detail.  Posthumous recognition was finally granted to Belousov with the award of the 1980 Lenin Prize.  This reaction, now known as the Belousov-Zhabotinskii or BZ reaction, is now considered the prototypical chemical oscillator and is still studied actively.

In this notebook we investigate limit-cycle oscillations using a simplified model due to Field and Noyes (1974) that faithfully describes the most important aspects of the relatively complicated chemical processes that occur in the BZ reaction.  The FN model is formulated in terms of a system of 3 differential equations

ε x'[t]==(1 - x[t]) x[t] - y[t] x[t] + q y[t]
δ y'[t]==-q y[t] - x[t] y[t] + 2 f z[t]
z'[t]==x[t] - z[t]

where x and y represent the concentrations of malonic acid (Subscript[HBrO, 2]) and bromide ions \(Br\^-\), while z represents the ratio between the concentrations of the (Ce^(4+)) and (Ce^(3+)) states of the cerium catalyst.  The parameters ε, δ, and q are related to reaction rates while the parameter f is stoichiometric and is generally taken as 0.5.  Periodic oscillations of the catalyst concentrations, coupled to appropriate dyes, produce dramatic oscillations of the color of the mixture.  Standing wave patterns can also be exhibited but are beyond the scope of this presentation.

References:

There is now a vast literature on the BZ reaction and related topics.  We base our discussion upon the pertinent chapter of Murray's textbook on Mathematical Biology and include a few references to original work below.  A larger set of references can be obtained from Murray or any recent review article.

B.P. Belousov (1951): A periodic reaction and its mechanism — English translation included in Oscillations and Travelling Waves in Chemical Systems, edited by R.J. Field and M. Burger (Wiley, New York, 1985).

B.P. Belousov (1959): An oscillating reaction and its mechanism, (Sborn. referat. radiat. med., Medgiz, Moscow), p. 145.

A.M. Zhabotinskii, Periodic processes of the oxidation of malonic acid in solution, Biofizika 9, 306 (1964).

A.M. Zhabotinskii and A.N. Zaikin, Autowave processes in a distributed chemical system, J. theor. Biol. 40, 45 (1973).

R.J. Field and R.M.  Noyes, Oscillations in chemical systems, J. Chem. Phys. 60, 1877 (1974).

J.D. Murray, Mathematical Biology (Springer-Verlag, Berlin, 1989), ch. 7.

Initialization

ClearAll["Global`*"];
Off[General::spell, General::spell1];
$TextStyle = {FontFamily -> "Times", 
          FontSize -> 12};

$DefaultFont={"Times", 12}; 
$TextStyle={FontFamily->"Times", FontSize->12,FontSlant->"Italic"};

Needs["Utilities`Notation`"];
Needs["Graphics`Graphics`"]

Off[Symbolize::"boxSymbolExists"];

Symbolize[(Subscript[x, s])]; Symbolize[(Subscript[y, s])]; Symbolize[(Subscript[z, s])];
Symbolize[(Subscript[x, 0])]; Symbolize[(Subscript[y, 0])]; Symbolize[(Subscript[z, 0])]

concentrations[t] = {x[t], y[t], z[t]};
rates[t] = D[concentrations[t], t];

FNrates[
          t]  =  {(ε^(-1)) (q y[t] - x[t] y[t] + 
                x[t] (1 - x[t])), (δ^(-1)) (-q y[t] - 
                x[t] y[t] + 2 f z[t]), x[t] - z[t]};

FNeqs  =  MapThread[Equal, {rates[t], FNrates[t]}]

TypicalParameters  =  {f -> 0.5, δ -> (1/3), 
          q -> 5×(10^(-3)),  ε -> 
            0.01};

Stability of steady states

Steady state

The following rules represent steady-state conditions in which the concentrations {x[t], y[t], z[t]} approach the steady values {(Subscript[x, s]), (Subscript[y, s]), (Subscript[z, s])}.

SteadyStateRules  = {x[t] -> (Subscript[x, s]), y[t] -> (Subscript[y, s]), 
          z[t] -> (Subscript[z, s]), (∂_t) x[t] -> 
            0, (∂_t) y[t] -> 0,  (∂_t) z[t] -> 
            0} ;

Applying these rules yields a set of equations

SteadyStateEquations = FNeqs  /. SteadyStateRules

with solutions

SteadyStates  = 
      Solve[SteadyStateEquations, {(Subscript[x, s]), (Subscript[y, s]), (Subscript[z, s])} ]  // 
        Simplify

SteadyStates /. TypicalParameters

The first is trivial while the second yields negative concentrations and is unphysical, leaving the third solution as the relevant case.

SteadyState = 
        SteadyStates[[3]];

Stability matrix

stabilityMatrix = 
      Outer[D, FNrates[t] /. SteadyStateRules, {(Subscript[x, s]), (Subscript[y, s]), (Subscript[z, s])}]

stabilityMatrix = stabilityMatrix /. SteadyState // Simplify

Although it is possible to evaluate the eigenvalues symbolically, the resulting expressions are rather unwieldy and the simplification process can exhaust the memory of typical desktop computers.  Therefore, we shall be content to give numerical values for typical parameters.

Eigenvalues[stabilityMatrix /. TypicalParameters]

Limit-cycle oscillations

The following function produces and displays solutions to the Field-Noyes model of the BZ reaction.

\(\(\(DisplayBZ::usage = \*""DisplayBZ[parameters,{(Subscript[x, 0]),(Subscript[y, 0]),(Subscript[z, 0])},{(Subscript[x, min]),(Subscript[y, min]),(Subscript[z, min])},{(Subscript[x, max]),(Subscript[y, max]),(Subscript[z, max])},(Subscript[t, max])] displays concentration variables for the Field-Noyes model of the BZ reaction for time interval {t,0,(Subscript[t, max])} assuming initial values {(Subscript[x, 0]),(Subscript[y, 0]),(Subscript[z, 0])}.  The list of parameters must give replacement rules for ε, δ, q, and f.  Minimum and maximum values are given for the plot ranges."";\)\(
\)
    \)	\)
\(DisplayBZ[
        parameters_, {\(x\_0\) : _, \(y\_0\) : _, \(z\_\(\(0\)\(\ \)\)\) : _}, {xmin_, ymin_, zmin_}, {xmax_, ymax_, zmax_}, tmax_] := 
      Module[{sol, plotx, ploty, plotz}, 
        sol = \(NDSolve[
              Join[FNeqs, {x[0] == x\_0, y[0] == y\_0, 
                    z[0] == z\_0}] /. parameters, {x, y, z}, {t, 0, 
                tmax}, MaxSteps -> 
                10\^4]\)[[1]]; 
plotx = 
          LogPlot[Evaluate[x[t] /. sol], \ {t, 0, tmax}, 
            PlotPoints -> 100, PlotRange -> {All, {xmin, xmax}}, \ 
            PlotLabel -> \*""(Subscript[HBrO, 2])"", 
            DisplayFunction -> Identity]; 
        ploty = LogPlot[Evaluate[y[t] /. sol], \ {t, 0, tmax}, 
            PlotPoints -> 100, PlotRange -> {All, {ymin, ymax}}, 
            PlotLabel -> \*""\!\(\(Br\^-\)\)"", 
            DisplayFunction -> Identity]; 
        plotz = LogPlot[Evaluate[z[t] /. sol], \ {t, 0, tmax}, 
            PlotPoints -> 100, PlotRange -> {All, {zmin, zmax}}, 
            PlotLabel -> \*""(Ce^(4+))/(Ce^(3+))"", DisplayFunction -> Identity]; 
        Show[GraphicsArray[{{plotx}, {ploty}, {plotz}}], Frame -> True, 
          FrameLabel -> {"t  (arb. units)", ""}, 
          PlotLabel -> "BZ reaction", 
          DisplayFunction -> $DisplayFunction]]\)

For example, even when the initial conditions are rather far from the unstable fixed point, the system settles upon a stable limit cycle after the sometimes dramatic transients have decayed.  The appearance of a periodic oscillation after an extended interval of monotonic variation is quite remarkable.

DisplayBZ[
        TypicalParameters, {0.0, 1.0, 1.0}, {0.001, 0.5, 0.001}, {0.5, 5.0, 
          0.5}, 100];