The law of mass action Consider a chemical reaction k A+B a C, where A, B and C represent chemicals, A and B are the reactants, C' is the product, and k > 0 is the rate constant of the reaction. The law ofmass action describes how the concentrations of A, B and C change as a consequence ofthis reaction. The idea is that the reaction will take place if chemicals A and B collide, and if their energy is higher than the energy of activation of the reaction. The number of successful (i.e. leading to a reaction) collisions between A and B is proportional to the product of the concentrations of A and B. The constant of proportionality is the constant k. We thus write the following differential equations for the expected values of [A] [B], [C], (re- call that we are in a situation where large numbers of molecules of reactants and prod- ucts are present) diA] diBl dici 7 : 7 : 'klAl [Bl : '7' where [1] denotes the concentration of chemical I. If A : B, so that the chemical re- action is k 2A~>C, we need two molecules ofAto create one product molecule 0', leading to 1 d[A] 2 d[Cl ,7 : 7k A : ,7. 2 dt [ ] dt If a reaction is reversible, for instance if 1'51 A + B : C, then A is consumed at relative rate k1 and produced at relative rate k2, so that i : *kllAl [B] + MC], and similarly for [B] and [C]. Finally, if a reaction process involves more than one chemical reaction, the rates of change of a chemical due to each of the reactions are added up in order to obtain the global rate of change of that chemical. For instance, assume that we have a system de- scribed by 1"31 AWB : 0, 1'32 k3 CWA ~> D. Then the concentration of A evolves according to d[A] 7F:7MN@+MM7MMML and the rate equations for B, Cand D are %:7mum+mm, %;:MNM7MM7MMML d[D] 77:mmw. Note that the net reaction associated with this process is 2A + B > D. The rate equations must be written for each of the steps involved in the reaction, and not on the basis of the net reaction. Chemical reactions can thus be described in terms of coupled nonlinear ordinan/ differ- ential equations, and the theon/ of dynamical systems therefore applies to their analysis. The Brusselator and Oregonator models Because rate equations describing chemical oscillations are nonlinear, it is possible to observe reactions which are oscillatory in time, in the same way as the Lotka-Volterra equations may describe oscillations in a predator-prey system. The Belousov- Zhabotinsky reaction is the classical example of an oscillatory chemical reaction. When Russian biochemist Boris P. Belousov reported his findings in 1951, his results were initially received with disbelief. It is only after A.M. Zhabotinsky ] reproduced and im- proved Belousov's experiments, that the existence of oscillatory reactions was finally accepted. There are two classical models for oscillatory reactions, called the 'Brusselator" (proposed by a group in Brussels ]) and the "Oregonator" (proposed by chemists at the University of Oregon 4). We briefly discuss each of them below. The Brusselator Consider the following hypothetical chemical process (Glansdorff & Prigogine, 1971). k1 A -> X, K2 2X +Y -+ 3X, K3 B + X - Y+D, KA X -+ E, where the rate constants ki are all equal to 1. The corresponding rate equations for [X] and [Y] form the Brusselator model, which reads [x]p dt = [A] + [X]? [Y] - [B][X] - [X], (8.1) d Y] dt [X]' + [B][X],where [ A] and [ B] are parameters. It has a unique fixed point P, given by [X] = [A] and [Y] = [B]/[A]. From a dimensional analysis point of view, these expressions may look strange, but we lost track of the dimensions when we set all of the reaction con- stants (which had different dimensions) to unity (see exercises). 6 5 4 [Y ] 3 N 2 3 4 [X] Figure 8.1. Phase plane of system (8.1), with [A]=1 and [B]=3, plotted with the software PPLANE. The Jacobian of (8.1) about the fixed point P is J(P) = [B] - 1 [A]2 - [B - [A]2and its determinant is equal to [A]2 > 0. The stability of the fixed point therefore de- pends easy t iectori on the sign of the trace of J(P), which is equal to T : [B] 7 1 7 [A]? It is 0 see that if [X] or m are sufficiently large, then dM/dpq : 71, and tra- es are almost straight lines with slope -1. However, as [Y] gets close to 0, these trajectories cannot leave the first quadrant, since d[l/dt > 0 if [Y] : 0 and [) > 0. Since at that time we also have d[)/di 0). The second attractor is a limit cycle, towards which all traiectories converge. Only one trajectory is plotted in this figure. The corresponding plots of [IQ and m as functions of time are shown in Figure 8.2. One can see that on the limit cycle, [X] remains small most of the time, then quickly increases to a maximum value and comes back towards zero, in a periodic fashion. As [X] increases, [Y] drops abruptly, and then slowly grows back to its maxi- mum value while [X] remains small. Such relaxation oscillationsare typical of many os- cillatory chemical reactions. The Oregonator The Oregonatorwas proposed by R.]. Fields and RM. Noyes in 1974. It corresponds to the following chemical reactions k1 AWY a X, 162 XWY a P, k3 B+X a 2X+z, 1'74 2X ~> Q, 1'75 z ~> fY, where [A] and [B] are constant, and f is a stoichiometric factor. The corresponding rate equations are d X] = ki [A] [Y] - K2 [X] [Y] + k3 [B] [X] - 2k4 [X]2, dt d[ Y] -ki [ A] [ Y] - k2 [ X] [Y] + ksf[Z], dt d[ Z] = k3 [B] [X] - ks [Z]. dt Field and Noyes (1974) defined the following dimensionless quantities [X] [Y] [Z] t C= Z= T = XO ' y = Yo Zo To where 1 Xo = K2 K1 [A], YO = K3 B, Zo = kik3 3 [ AB], To = K 2 K2ks VKiks [A][B] Then, the dimensionless form of the Oregonator becomes dx dT =s(y - xy + x - qx2), dy dT = =(-y -xy+ fz), (8.2) dz dT = w(x - z), where K3 B 2kik4[A] k5 S = q= W = KI[A]' K2k3[B] VKik3 [A][B] By assuming that a quickly reaches a steady-state value, it is possible to reduce Equations (8.2) to a two-dimensional dynamical system. Indeed, setting da /dT = 0 givesC = 2q - y+ ( 1 - y)2 + 4qy) and by substituting this expression into (8.2), one obtains dy = ( -y - y dT 2q (1 - y+ (1 - y)2 + 4qy ) + fz), (8.3) dz = W dT 1 - y + ( 1 - y) 2 + 4qy) - z) . 35 30 25 20 15 10 5 2 4 5 6 w y Phase plane of system (8.3), with f=1, 5=1, q=0.01 and w=2, plotted with the software PPLANE.In what follows, we set f : 1. Field and Noyes (1974) estimated the values of s, q, u) as well as of all dimensionless variables for the BelousovZhabotinsky reaction. Details can be found in their article entitled Oscillations in chemical systems. IV Limit cycle be- havior in a model of a real chemical reaction. They also concluded that the Oregonator model is capable of predicting oscillations for the concentrations of the chemicals in- volved in the reaction process. Realistic parameters make the problem very stiff (i.e. there are both ven/ short and very long characteristic time scales, which is typical for excitable systems), but model () also exhibits oscillations in reasonably stiff cases. Figure 8.3 shows the phase plane of model (g) with parameters 5 : 1, q : 0.01, and w : 2. Trajectories starting from the fixed point with both 31 and z non-zero con- verge towards a limit cycle. As a consequence, the variables y and z oscillate in time, and so do [)q, [Y] and [Z] in the Oregonator model. Chemical waves When a chemical reaction takes place in a spatially extended system, diffusion should be taken into account. As a consequence, the rate equations discussed above are turned into partial differential equations. The reaction terms remain the same, but the concentration of each chemical now diffuses with a diffusion coefficient which de- pends on the size and weight of the molecule in question. If the reaction is oscillatory, wave fronts, corresponding to say large concentrations of some chemical, propagate in the system. Moreover, if the reaction is constrained to a twodimensional surface and the wave is initiated at some point in space, then the wave fronts are circular. As an ex- ample, a 2001 article by C. Sachs et ail5| describes how such wave fronts are observed in an experiment, and proposes a reaction-diffusion model which reproduces this be- havior. The authors also discuss the limitation of reactiondiffusion models when macroscopic parameters are affected by the details of microscopic interactions. What would happen if the wave front was broken, for instance if the wave went over a region where the reaction could not take place? If the wave front is anchored at one point, then it will curve and eventually form a spiral wave. Such waves are often ob- served in experiments where the Belousov-Zhabotinsky reaction is constrained to a two-dimensional surface, such as a thin film, a porous glass disk, or even a sheet of fil- ter paper. The article by SC. Miller et allEl discusses experiments revealing the struc- ture of the core of such spiral waves. Summary We started with the law of mass action, which allowed us to describe the dynamics of the (average) concentrations of reactants and products involved in chemical reactions. In particular, we considered 'two hypothetical sets of chemical reactions called the Brusselator and the Oregonator. By applying the methods of analysis discussed in the previous chapters, we concluded that it was possible for these chemical systems to ex- hibit oscillatory dynamics. This provides a proof of concept for oscillatory reactions such as the Belousov-Zhabotinsky reaction. In spatially extended systems, the latter leads to chemical waves and spiral defects, whose structure is well documented in the research literature. Problem 3 Consider the chemical reaction k nA -+ B, where n is a positive integer. 1. What is the rate equation for A? 2. Given that [A] (0) = [A]o, find the solution of this equation. 3. Does the solution make sense for all times? If not, what happens