J. Phys. Chem. 1996, 100, 8323-8330
8323
The Use of Delay Differential Equations in Chemical Kinetics Marc R. Roussel* Centre for Nonlinear Dynamics, McGill UniVersity, Montreal, Canada ReceiVed: January 4, 1996X
The concept of a chemically acceptable model is developed. Chemically acceptable models are causal and maintain the nonnegativity of concentrations. An extension of the law of mass action allowing delayed effects is described and shown to lead to chemically acceptable models. Delayed variable enzyme catalysis and Oregonator models are studied and shown to be equivalent to their classical mass-action counterparts from a number of perspectives. This equivalence is shown to be a consequence of an exact correspondence between certain ordinary and distributed-delay differential equations. For sufficiently small delays, this correspondence extends to discrete-delay systems. However, reversible chemical reaction networks cannot be adequately modeled with discrete-delay equations.
1. Introduction In several recent chemical modeling studies, delay differential equations (DDEs) have been used instead of the usual massaction ordinary differential equations (ODEs).1-6 There are well-established laws regarding the class of ODEs which can represent a chemical system,7-10 but little has been written about the possibility of using other classes of equations. Delay differential equations present both significant opportunities and unique difficulties for chemical modeling. There are many reasons why we might wish to expand the class of models available to chemists to include delayed variable formulations. The emphasis in creating a DDE model of a chemical system is shifted from cataloging intermediates and their reactions to describing the dynamic relationships between the concentrations of key species. As a result, fewer concentration variables will generally appear than in a classical massaction mechanism. In preliminary investigations of a new mechanism, the level of description afforded by a DDE model is often closer to our state of knowledge than is a detailed mechanism in which a certain amount of speculation about intermediate species is a necessary element. Despite these advantages, some caution should be exercised in using delayed variable models in chemistry. The evolution of closed chemical systems is subject to nontrivial constraints which in turn have strong implications for the form of realistic models of open systems. Some of these issues are examined in this paper. Examples of successful DDE models are presented, but some of the limitations of such models are also assessed. A context is provided for this discussion by the introduction of a set of axioms which a very general class of chemical models must satisfy (section 2) and of a delayed variable extension to the law of mass action (section 3). A detailed study of Brown’s model of enzyme catalysis follows. The computation of the equilibrium point of a DDE model is not entirely straightforward. The procedure is discussed, using Brown’s model as an example, in section 4. Chemical relaxation usually has a hierarchical nature with faster and slower processes,11-13 leading to the formation of slow invariant manifolds.14-16 In section 5 we discover such a * Permanent address: Department of Chemistry, University of Lethbridge, Lethbridge, Alberta, Canada T1K 3M4. E-mail: roussel@ henri.chem.uleth.ca X Abstract published in AdVance ACS Abstracts, May 1, 1996.
S0022-3654(96)00067-6 CCC: $12.00
structure in the solution space of the Brown model. Oscillations around the slow manifold are observed, as they are in certain related ODE models. In a recent paper, Epstein and Luo studied two delayed variable Oregonators.3 Building on the delayed mass-action framework introduced in section 3, a third version is presented and analyzed in section 6. This delayed Oregonator reproduces the waveforms of the ODE model surprisingly accurately for small delays. The close correspondence between the Oregonator and delayed Oregonator is rationalized in section 7 by showing that the Oregonator is equivalent to a distributed-delay system of form similar to the delayed Oregonator studied here. The argument is quite general and could be used to eliminate linear subsystems in any mechanism. Finally, in section 8, two generalizations of Brown’s model are studied. The first example shows the kinetic equivalence between Brown’s model and a similar system in which the enzyme does not immediately return to a reactive conformation on release of the product. In the second example, it is shown that naive attempts to incorporate reversibility into Brown’s model lead to behavior in violation of the law of microscopic reversibility. This implies that simple delay systems cannot in general be used to model chemical behavior near the equilibrium point. 2. Chemically Acceptable Models The law of mass action generates ODE problems of the initialvalue type; that is, in addition to the differential equations themselves, the concentrations at t ) 0 are required. In DDE models, it is necessary to specify an initial function for each of the dependent variables in the range t [-τmax, 0] where τmax is the largest of the delays in the system.17 In the case of isothermal chemical kinetics models, the dependent variables are the concentrations. To guarantee that the solutions of the DDEs remain nonnegative, it will be necessary to constrain both the starting functions and the class of allowable DDEs.18 In ODE models, this is accomplished by invoking some form of the law of mass action. The classical law of mass action8 is a little restrictive and excludes important classes of chemical models.10 Let us define a chemically acceptable model: Let x(t) Rn be the vector of chemical concentrations in the system under study at time t; the elements of x are denoted xi. The evolution equation © 1996 American Chemical Society
8324 J. Phys. Chem., Vol. 100, No. 20, 1996
Roussel products, the rate of formation of any species is never negative on the boundary of the physically realizable part of concentration space.
is of the form
x3 (t) ) f(x(τ),t) where f is a vector-valued functional with elements fi. The full model includes an initial vector-valued function O(t) (elements φi(t)) defined in the interval t [-τmax, 0], where τmax is the greatest delay in the system.19 A model shall be said to be chemically acceptable if (1) f depends on x only up to time t; (2) fi g 0 whenever xi(t) ) 0; (3) φi(t) g 0 for t [-τmax, 0]. Rule 1 disallows terms which depend on the future values of the variables and thus ensures that causality is not violated. Rules 2 and 3 guarantee that negative concentrations will not be generated by the model. The class of chemically acceptable models includes as special cases all previous formulations of generalized mass-action laws.7-10 It is important to note that these are minimal constraints. Important properties of chemical systems such as microscopic reversibility are not captured by these rules. This is partly deliberate: Far-from-equilibrium chemical models frequently violate near-equilibrium constraints. This is acceptable as long as we keep in mind that such models are only valid in certain dynamical regimes. These questions will be examined in greater detail in section 8.2.
In 1902, Brown proposed a model of enzyme catalysis in which the enzyme-substrate complex, rather than decaying by the usual first-order rate law, has a fixed lifetime.20 This model of enzyme kinetics was quickly succeeded by Henri’s massaction treatment.21 It is nevertheless interesting to investigate the kinetics of Brown’s model and to inquire whether any realistic mass-action systems are (in some sense) equivalent to it. Before we proceed, it is convenient to devise a new notation for systems such as that described by Brown. In summary, Brown describes some ordinary chemistry (complex formation) followed by a delay, concluding with the appearance of some products. Such delayed mass-action steps are denoted by
where k is a rate constant, {τi} is the set of delays (each of which is nonnegative), and {Ri} and {Pi} are the sets of reactants and products, respectively. When a step of this sort is written, it will be understood that the disappearance of the reactants follows the usual mass-action law. Each product Pi, on the other hand, appears τi time units after the mass-action step. Brown’s model, stated in this notation, is exceedingly simple:
(Brown)
The DDEs generated by this model are
(2)
(This equation is a simple consequence of the stoichiometry of the reaction and is equally true for the usual Michaelis-Menten mechanism or indeed for any other network of the form
E + S h ... h P + E regardless of the contents of the ellipsis.) The DDEs 1 are most easily handled in dimensionless form. Accordingly, let
θ ) t/τ,
y(θ) ) kτe(t)
x(θ) ) kτs(t), z(θ) ) kτp(t) The Brown system then becomes
y˘ (θ) ) -x(θ) y(θ) + x(θ-1) y(θ-1)
(3)
z(θ) ) µ - x(θ) + y(θ) where dots now represent differentiation with respect to the scaled time θ and the redundant DDE has been replaced using eq 2, defining
µ ) kτ∆m By definition, equilibrium solutions satisfy x(θ) ) x(θ-1) and y(θ) ) y(θ-1), which makes y˘ identically zero. The equilibrium condition therefore reduces to xeqyeq ) 0, but since the mechanism is irreversible, the solution is obvious:
In an ODE model, it would be necessary to refer to the initial concentrations to determine yeq and zeq. Similarly, we shall now discover that the coordinates of the equilibrium point depend on the initial functions in a DDE model. From the conservation relation, we know that zeq - yeq ) µ, so that we need only find one to know the other. Consider the following gedanken experiment: Imagine that at time θ0 the reaction is quenched (e.g., the enzyme is poisoned), rapidly producing an inactive enzyme (scaled concentration jy) from any free enzyme present in the system. Then, y obeys
jy˘ (θ) ) x(θ-1) y(θ-1)
s˘ (t) ) -ke(t) s(t) e˘ (t) ) -ke(t) s(t) + ke(t-τ) s(t-τ)
s + p - e ) ∆m
xeq ) 0
k
R1(t) + R2(t) + ... 98 P1(t+τ1) + P2(t+τ2) + ...
k
First note that one of the equations of set 1 is redundant since s˘ + p˘ - e˘ ) 0, so that, at any time,
x˘ (θ) ) -x(θ) y(θ)
3. Delayed Mass-Action Systems
E(t) + S(t) 98 P(t+τ) + E(t+τ)
4. The Equilibrium Point of Brown’s Model
(1)
p˘ (t) ) ke(t-τ) s(t-τ) (As usual, overdots represent differentiation with respect to time.) It is easy to verify that systems consisting uniquely of ordinary and delayed mass-action steps are always chemically acceptable (as defined in the previous section), provided appropriate initial functions are constructed: Advanced terms are not generated, and since delayed terms only represent the appearance of
for θ [θ0,∞] with jy(θ0) ) y(θ0). After θ ) θ0 + 1, all reaction stops and
jy(θgθ0+1) ) y(θ0) + ∫θ
θ0+1
x(ζ-1) y(ζ-1) dζ
0
The inactive form eventually collects all the available enzyme so we can guess that this is just the total (scaled) amount of enzyme, a quantity that should be conserved. In fact, it is easy to verify that
ytotal ) y(θ) + ∫θ-1x(ζ) y(ζ) dζ θ
Use of Delay Differential Equations in Chemical Kinetics
J. Phys. Chem., Vol. 100, No. 20, 1996 8325
is a constant of the motion. In particular, evaluating this constant at equilibrium and at θ ) 0, we have
yeq ) y(0) + ∫-1x(ζ) y(ζ) dζ 0
If z(0) ) 0, the equilibrium point is therefore
(xeq, yeq, zeq) ) (0, y(0) + M, z(0) + M)
(4)
where
M ) ∫-1x(θ) y(θ) dθ 0
5. Global Analysis of Brown’s Model If systems of DDEs are to be used as exploratory tools in chemical modeling, they must mimic the phase space structure of ODE models. Brown’s model is a DDE counterpart of the classical Michaelis-Menten mechanism, whose phase space structure is well-known: Trajectories starting from any point in the physically realizable part of phase space are attracted to the equilibrium point along a unique curve,22 the slow manifold M.23-25 Generalizations of the Michaelis-Menten mechanism with multiple intermediate complexes can display transient oscillations about M,26-28 but almost every physically realizable trajectory approaches the equilibrium point along this manifold.23,29,30 In this section, the Brown model is shown to possess an analogous property: For almost every positive initial function with a sufficiently large value of x(0), the system executes transient oscillations about the slow manifold; the equilibrium point is globally attractive and is approached along M. Although we could immediately launch into a mathematical analysis, it is easiest to begin by examining trajectories of the model equations obtained by numerical integrations.31 There is an important point to note here: Although trajectories are displayed in an x × y “phase plane”, we cannot think of this plane as we do in ODE models because incomplete information is being presented in the DDE case, the dependence on history being, in a sense, hidden in these pictures. The behavior of the model with a variety of initial functions was examined. Two different normalization procedures for the initial functions can be used, corresponding to two different classes of experiments.32 1. Values of yeq and zeq are chosen and the initial functions are normalized in such a way that the equilibrium relations 4 are obeyed. This corresponds to experiments with a fixed equilibrium composition. 2. The y component of the initial function is normalized in such a way that yeq is the same in a given series of numerical experiments, but the total amount of substrate is left as a free parameter. This corresponds to experiments in which a fixed amount of enzyme processes different amounts of substrate; this type of study is most appropriate when the mechanism is irreversible or nearly so. The second normalization procedure is the more revealing in this case. Figure 1 shows several trajectories for the Brown model started from random initial functions,33 along with the slow manifold, whose numerical construction is discussed below. As long as we start from sufficiently far out along the x axis, trajectories started from almost any initial function relax to M in an oscillatory fashion. In fact, of all the initial functions tried, only ones which are very nearly constant at large x did not produce oscillations. The amplitude of the oscillations depends on the starting function, but the frequency is a characteristic property of the DDEs.
Figure 1. Several trajectories for the Brown model. The dotted trajectory is the numerically computed slow manifold. The others were started from random initial functions. After a short time, all but one oscillatory mode is filtered out, which then decays, leaving the system to approach the equilibrium point at (0, 1) along the slow manifold.
The slow manifold is constructed numerically by starting a trajectory from a constant initial function, normalized to achieve the correct equilibrium value of y, with a sufficiently large value of x(0).34 In other words, a representation of the slow manifold valid at large x can be recovered from the normalization condition for a constant initial function (eq 4, specialized to the case of constant functions):
yeq ) yM + xyM or
yM(x) ≈
yeq 1+x
(5)
This is exactly the formula for the steady-state approximation in the Michaelis-Menten mechanism, which is known to be an accurate formula for the slow manifold of that mechanism at large substrate concentrations!24 We can also compute the slope of the slow manifold at the equilibrium point by linearization of eqs 3: Let
δx(θ) ) x(θ) - xeq ) x(θ) and
δy(θ) ) y(θ) - yeq Then, in the vicinity of the equilibrium point
δx˘ (θ) ) -yeqδx(θ), δy˘ (θ) ) yeq(δx(θ-1) - δx(θ)) The first of these equations has a trivial solution:
δx(θ) ) x(θ0)e-yeq(θ-θ0) Substituting this solution into the second differential equation and solving, we get
δy(θ) ) x(θ0)e-yeq(θ-θ0)(1 - eyeq) Thus, all solutions approach the equilibrium point with slope
δy/δx ) 1 - eyeq
8326 J. Phys. Chem., Vol. 100, No. 20, 1996
Roussel
and, in particular, this must be the slope of the slow manifold there. Note that since yeq > 0, the slope is negative, as it must be. It is remarkable not only that this model possesses a slow manifold but that this manifold has a geometry similar (as witnessed by eq 5) to that of the Michaelis-Menten mechanism. Some reasons for this close equivalence will be discussed in a later section of this paper. Note that it may be possible to adapt the methods of Wischert et al.35 to the approximation of invariant manifolds in delayed mass-action systems. 6. A Delayed Oregonator Model The classical Oregonator model of the Belousov-Zhabotinskii reaction is36,37 Figure 2. Numerical solutions of the Oregonator (s) and delayed Oregonator (‚‚‚) equations for the FKN parameter set36 R ) 77, β ) 0.16, f ) 1.0, γ ) 12, ) 0.025. In this and subsequent figures, the Oregonator initial condition is (x, y, z) ) (1, 1, 1); the DO system initial function is (x(θ), y(θ)) ) (1, 1) for θ [-, 0]. The two solutions are an excellent qualitative match.
k1
A + Y 98 X + P k2
X + Y 98 2P k3
A + X 98 2X + Z
(Oregonator)
k4
2X 98 A + P k5
Z 98 fY It has been claimed that the only role of Z is to introduce delayed feedback.3,38 This seems reasonable given the manner in which Z appears in the network and suggests replacing the third and fifth steps of the Oregonator by the delayed mass-action step k3
A(t) + X(t) 98 2X(t) + fY(t+τ) We thus arrive at the following delayed Oregonator (DO) model: k1
A(t) + Y(t) 98 X(t) + P(t) k2
X(t) + Y(t) 98 2P(t)
(DO)
k3
A(t) + X(t) 98 2X(t) + fY(t+τ)
x˘ ) R[fβy - fRxy + Rx - 2βx2]
while the DO model obeys the equations
x˘ (θ) ) R[fβy(θ) - fRx(θ) y(θ) + Rx(θ) - 2βx(θ)2] y˘ (θ) ) -y(θ) - (R/β)x(θ) y(θ) + (R/β)x(θ-)
k5 k1A
θ ) k1At
γ)
) k1Aτ
x(θ) )
X A
R ) xk3/k1
y(θ) )
Y k2 A fk3
x
Z k5 A k3A
(6)
z˘ ) γ(x - z)
The concentration of A is taken to be fixed by replenishment from a large external reservoir; P is a sink species whose kinetics we can ignore. We proceed immediately to the dimensionless equations for the two systems by defining39
z(θ) )
The Oregonator ODEs are then
y˘ ) -y - (R/β)xy + (R/β)z
k4
2X(t) 98 A(t) + P(t)
β ) xk4/k2
Figure 3. Numerical solutions of the Oregonator (s) and delayed Oregonator (‚‚‚) equations for a set of parameters inspired by the FieldFo¨rsterling set37 with a moderate value of γ but f somewhat smaller than normal; R ) 5.1, β ) 0.035, f ) 0.57, γ ) 2.0, ) 0.84. The x components of the solutions agree reasonably well; the y components less so.
k2k4 k1k3
x
k2k4 k1k3
(7)
The rate constants for the Oregonator are known.37 What of the delay τ appearing in the DO model? Roughly speaking, τ is the mean time required for 1 mol of Z to be converted to f mol of Y. Thus, one obvious estimate for τ might be the halflife of Z, i.e. τ ≈ ln 2/k5 or, in dimensionless units, ≈ ln 2/γ. Epstein and Luo have proposed that in such systems an estimate of the delay can be computed as the lag between a maximum in the eliminated variable (z) and the next maximum in the delayed variable (y).3 Figures 2-4 show numerical simulations40 of the Oregonator and DO models in which was chosen to match the frequency of the Oregonator oscillator. The waveform is relatively
Use of Delay Differential Equations in Chemical Kinetics
Figure 4. Numerical solutions of the Oregonator (s) and delayed Oregonator (‚‚‚) equations for the same parameters as Figure 3 except for a much smaller value of γ and, correspondingly, a larger value of ; γ ) 0.26, ) 13.2. The larger the delay, the less well the solutions of the two models agree.
insensitive to . Also note that the initial function in the DO model has no effect on the eventual kinetics. When the two models produce qualitatively similar kinetics (Figures 2 and 3), either the half-life or Epstein and Luo’s lag time criterion can be used to provide an order-of-magnitude estimate of . For instance, with the Field-Ko¨ro¨s-Noyes (FKN) parameter set,36 the optimal delay is ) 0.025; from the half-life we compute ≈ 0.058, while the Epstein-Luo estimator gives ≈ 0.18. With the Field-Fo¨rsterling (FF) parameters used to produce Figure 3, the Epstein-Luo estimator is ≈ 0.82, while the optimal value is ≈ 0.84; the half-life estimator is ≈ 0.35. A third method for estimating relies on linearization. The two models have identical fixed points, one of which is an unstable steady state at the origin. The only other physically realizable fixed point is (x*, y*, z*):
x* ) z* ) βy*/[R(1 - y*)]
J. Phys. Chem., Vol. 100, No. 20, 1996 8327 {1440, 0.1010, 0.07521 ( 75.87i, -0.001787 ( 151.7i, ...}. Since the two leading eigenvalues match to within a fraction of a percent, this method also yields an order-of-magnitude estimate of . This case was relatively easy since the two leading eigenvalues are both real. When the leading eigenvalues are a complex-conjugate pair, the correct procedure is to match the imaginary parts (the frequencies). Quite apart from their use in estimating the size of , the above calculations show a part of the reason that the two models produce similar behavior: Their leading eigenvalues are similar. The importance of this property has previously been emphasized by Busenberg and Travis.41 Further notes on the equivalence between the two models are presented in the next section. Figure 4 shows a case in which the solutions of the Oregonator and DO models are really quite different. In this case, γ is relatively small, which, due to the inverse proportionality between γ and , implies that the delay is large. In general, we find that the two models agree best when the delay is small. 7. Notes on the Relationship between Classical and Delayed Mass-Action Systems Both the Brown and DO models were obtained from their classical counterparts by eliminating a single intermediate in favor of a delayed term. To better understand the relationship between classical and delayed mass-action systems, at least in these simple cases, we consider an exact reduction of the Oregonator to a distributed delay system. An identical argument could be used to transform the Michaelis-Menten mechanism into a distributed delay system. The third of the Oregonator ODEs 6 has the solution
z(θ) ) ∫0 x(ξ)γe-γ(θ-ξ) dξ + z(0)e-γθ θ
R2(1 - f) - 2β2 + x[R2(1 - f) - 2β2] + 8R2β2(1 + f)
After a sufficiently long time, the trailing term becomes negligibly small, provided we are not interested in transient behavior. Thus, the Oregonator ODEs are equivalent to
R2(1 - f) + 2β2 + x[R2(1 - f) - 2β2] + 8R2β2(1 + f)
x˘ (θ) ) R[fβy(θ) - fRx(θ) y(θ) + Rx(θ) - 2βx(θ)2]
y* ) 2
2
The characteristic equation of the Oregonator at this point is
βΛ3 + Λ2[Rβ(fRy* + 4βx* - R) + β + Rx* + βγ] + Λ[R(fRy* + 4βx* - R)(β + Rx* + βγ) + γ(β + Rx*) + fR2y*(β - Rx*)] + Rγ[(fRy* + 4βx* - R)(β + Rx*) fR(β - Rx*)(1 - y*)] ) 0 (8) The characteristic equation of the delayed Oregonator at the nontrivial fixed point is
y˘ (θ) ) -y(θ) - (R/β)x(θ) y(θ) + (R/β)∫0 x(ξ)γe-γ(θ-ξ) dξ θ
These equations are of a form similar to the DO model eqs 7 except that the delay is distributed. Again, discarding information about the transient behavior of the original model, for sufficiently large θ, the second of these equations is obviously equivalent to
y˘ (θ) ) -y(θ) - (R/β)x(θ) y(θ) + (R/β)∫-∞x(ξ)γe-γ(θ-ξ) dξ 0
βλ2 + λ[Rβ(fRy* + 4βx* - R) + β + Rx*] + R[(fRy* + 4βx* - R)(β + Rx*) + fRy*(β - Rx*)] ) fR2(β - Rx*)e-λ (9) An estimate of can then be obtained by matching up the leading eigenvalues of the two models. In cases in which the two models display qualitatively similar behavior, this method is about as accurate as the other two. With the FKN parameter set, we get Λ ) {1433, 0.1010, -1964}. At the value of which reproduces the frequency of circulation around the limit cycle, namely, ) 0.025, the leading eigenvalues of the DO model at the nontrivial fixed point are λ ) {1440, 0.3370, -0.6095 ( 253.1i, -3.304 ( 506.0i, ...}. For ) 0.083, λ )
(The conditions under which this extension of the range of integration is valid have been studied by Cooke and Grossman.42) Note that the delay distribution g1(u) ) γe-γu is the lowest order normalized gamma distribution.43 The normalization (∫∞0 g1(u) du ) 1) simply expresses mass conservation through the delay element. We can now see why the Oregonator and DO models agree for large γ (small ) but not otherwise: We can rigorously write the Oregonator equations as a pair of distributed-delay equations. For small γ, the delay distribution is very broad and a single finite delay is simply not a good model. On the other hand,
8328 J. Phys. Chem., Vol. 100, No. 20, 1996
Roussel
for large γ, only a very small part of the past history is relevant so that a single delay can capture the essential features of the dynamics. Consider an irreversible linear chain (f C1 f C2 f ... f Cq f) in which all the rate constants are the same. Such a chain can be replaced by a term in the governing differential equations with a gamma-distributed delay.43 The delay distribution then has the form gq(u) ) kquq-1e-ku/(q - 1)!, where k is the common rate constant. As q increases, particularly for large values of k, gq becomes more sharply peaked. Thus, the Brown and DO models, with their single delays, are most closely equivalent to mechanisms with multiple intermediates, each of which is converted irreversibly to the succeeding form. Furthermore, it is known that multiple-complex Michaelis-Menten mechanisms can produce transient oscillations about their slow manifold and that the damping of these oscillations is minimized by making the reactions irreversible and setting all the forward rate constants in the chain to a common value.28 By the arguments presented above, the Brown model mimics a Michaelis-Menten network constrained in precisely this manner and the observation of oscillations about the slow manifold in the former is therefore not surprising. 8. Generalizations of Brown’s Model Additional insight into the utility and limitations of delayed mass-action systems may be gained by studying two generalizations of Brown’s model. 8.1. The Brown Model with Enzyme Recycling. Suppose that, after releasing the product, the enzyme finds itself in an unreactive conformation. Some time must then elapse for it to return to its native state. If we were to codify this hypothesis with a conventional ODE model, we would need an extra concentration variable (for the inactive enzyme) and an extra parameter (the rate constant for enzyme recycling). In the framework developed here however, we need only add a parameter to the Brown model to cover this case, which we call the Brown with recycle (BWR) model: k
E(t) + S(t) 98 P(t + τ1) + E(t + τ2)
(BWR)
equations of the Brown model (eqs 3). It follows that the BWR model is kinetically equivalent to the Brown system. This in turn suggests that enzyme recycling is unlikely to produce interesting dynamics when the basic reaction scheme is of the Michaelis-Menten type. The BWR model is a delayed variable version of a “Michaelis-Menten with recycle” (MMWR) mechanism,
E+ShCfF+P FfE which has equations of motion identical to an irreversible twocomplex Michaelis-Menten mechanism, namely,
E + S h C1 f C2 f E + P Thus, kinetic equivalence is preserved in the transformation of classical to delayed mass-action models. 8.2. The Reversible Brown Model. Enzymes are ordinary molecules whose reactions are subject to the law of microscopic reversibility.44 In modeling studies, it is often highly desirable to incorporate all the reverse steps of a mechanism.45 We might naively try to incorporate reversibility in the Brown model by adding a reverse reaction: k1
E(t) + S(t) 98 P(t+τ1) + E(t+τ1) k-1
E(t) + P(t) 98 S(t+τ-1) + E(t+τ-1)
A little reflection, however, reveals that this model is more closely related to
x(θ) ) kτ2s(t)
τ ) τ1/τ2
y(θ) ) kτ2e(t)
C1
C2
E+P
k1
k-2
-1
2
E + S y\ z C y\ zP+E k k
(RMM)
Since model RB1 is not microscopically reversible, its solutions oscillate about the equilibrium point. Classical mass-action models never display such behavior.46 A better approach starts from the RMM mechanism. We formally eliminate c by solving its differential equation and apply the Cooke-Grossman range extension42 to transform the resulting integro-differential equations into distributed delay equations, as was done in section 7 for the Oregonator. This leaves us with the equations
z(θ) ) kτ2p(t) from which it follows that
x˘ (θ) ) -x(θ) y(θ) y˘ (θ) ) -x(θ) y(θ) + x(θ-1) y(θ-1)
E+S
than to the desired reversible Michaelis-Menten (RMM) mechanism,
It is understood that τ2 g τ1. We introduce a scaling similar to that used in the Brown model:
θ ) t/τ2
(RB1)
(10)
z˘ (θ) ) x(θ-τ) y(θ-τ) This system obeys the unconventional conservation relation
c(t) ≈
∫-∞t e-(k
-1+k-2)ξ
[k1e(t-ξ) s(t-ξ) + k2e(t-ξ) p(t-ξ)] dξ
e˘ (t) ) -k1e(t) s(t) - k2e(t) p(t) + (k-1 + k-2)c(t)
z(θ) ) µ + y(θ+1-τ) - x(θ+1-τ)
s˘ (t) ) -k1e(t) s(t) + k-1c(t)
Since τ < 1, if we rely on this formula to compute z, our knowledge of the product concentration lags our data on x and y. Nevertheless, the first two equations of the set 10 are clearly sufficient to construct the kinetics of the full system. However, note that these equations are identical to the corresponding
p˘ (t) ) -k2e(t) p(t) + k-2c(t)
(11)
These equations should have solutions essentially identical to those of the reversible Michaelis-Menten mechanism, except perhaps for the early (transient) evolution. Next, we attempt
Use of Delay Differential Equations in Chemical Kinetics to replace the integral by a simple approximation. The kernel is exponentially decaying with a characteristic decay time (k-1 + k-2)-1. Therefore, for some mean delay τ,
c(t) ≈ [k1e(t-τ) s(t-τ) + k2e(t-τ) p(t-τ)]/(k-1 + k-2) We also note that s˘ + p˘ - e˘ ) 0. This leaves us with a pair of coupled delay differential equations, the RB2 model. Despite the more sophisticated balancing, the solutions of these equations also oscillate about the equilibrium point, which, as mentioned above, is forbidden by the law of microscopic reversibility.46 Our failure to generate a delayed mass-action model equivalent to the reversible Michaelis-Menten mechanism near the equilibrium point is a result of a general feature of the characteristic equations of delay differential equations: These characteristic equations generally have an infinite number of complex roots.17 As a result, we can expect damped oscillations about the equilibrium point of almost every reversible model incorporating simple delays. Chemically acceptable models of an irreversible reaction network, on the other hand, are unlikely to oscillate near the equilibrium point because this point will appear on the boundary of the physically realizable region and, by construction, their solutions cannot cross this boundary. The obvious solution to this problem is to halt our reduction process at the distributed-delay stage, i.e. to accept sets of equations such as 11. These equations satisfy all of the requirements of a chemically acceptable model, and they behave correctly near the equilibrium point since they are fully equivalent to the mass-action ODEs. 9. Discussion It appears from this study that simple delayed mass-action systems, used with caution, can be of some use in chemical modeling. The Brown and DO models, which replicate the behavior of well-studied mass-action models within their respective domains of validity, are witnesses to this. The most serious outstanding problem is the treatment of fully reversible mechanisms operating near equilibrium or near stable steady states satisfying detailed balance.10,47 These systems require properly distributed delays. The Brown system and its variations are models of closed chemical systems. A closed chemical system must have a unique, globally stable equilibrium point.47,48 Proving uniqueness of the equilibrium for delayed mass-action models representing closed systems is probably straightforward since the equilibrium is, by definition, homogeneous in time. The equilibrium relations therefore reduce to those of ordinary massaction systems, give or take the computation of the total amount of each species, as discussed in section 4. Shear’s proof of uniqueness48 can thus likely be extended to delayed mass-action systems. On the other hand, stability is likely to be a much more difficult problem. The classical proof of the linear stability of chemical equilibria46 probably has no straightforward extension to systems with delays because the type of problem changes from showing that a certain class of polynomials has no roots with positive real parts to showing the same thing but for a class of transcendental equations. Interestingly, extending Shear’s proof of global stability47 to the delayed mass-action case may be easier, although the application of Lyapunov function methods to systems of delay differential equations presents a number of technical difficulties.49 Acknowledgment. I would like to thank Prof. Simon Fraser for encouraging me to pursue this work. I would also like to thank Prof. Michael Mackey and Dr. Je´roˆme Losson for valuable
J. Phys. Chem., Vol. 100, No. 20, 1996 8329 discussions on the properties and analysis of DDEs. During this research, I was supported by a Postdoctoral Fellowship from the Natural Sciences and Engineering Research Council of Canada. References and Notes (1) Schell, M.; Ross, J. J. Chem. Phys. 1986, 85, 6489-6503. (2) Epstein, I. R. J. Chem. Phys. 1990, 92, 1702-1712. (3) Epstein, I. R.; Luo, Y. J. Chem. Phys. 1991, 95, 244-254. (4) Cheng, Z.; Hamilton, I. P.; Teitelbaum, H. J. Phys. Chem. 1991, 95, 4929-4931. (5) Chevalier, T.; Freund, A.; Ross, J. J. Chem. Phys. 1991, 95, 308316. (6) Bar-Eli, K.; Noyes, R. M. J. Phys. Chem. 1992, 96, 7664-7670. (7) Wei, J. J. Chem. Phys. 1962, 36, 1578-1584. (8) Horn, F.; Jackson, R. Arch. Rational Mech. Anal. 1972, 47, 81116. (9) Perelson, A. S.; Wallwork, D. J. Chem. Phys. 1977, 66, 43904394. (10) Schuster, S.; Schuster, R. Z. Phys. Chem. 1990, 271, 337-345. (11) Curtiss, C. F.; Hirschfelder, J. O. Proc. Natl. Acad. Sci. U.S.A. 1952, 38, 235-243. (12) Heinrich, R.; Rapoport, S. M.; Rapoport, T. A. Prog. Biophys. Mol. Biol. 1977, 32, 1-82. (13) van Kampen, N. G. Phys. Rep. 1985, 124, 69-160. (14) Okuda, M. Prog. Theor. Phys. 1981, 66, 90-100. (15) Fraser, S. J. J. Chem. Phys. 1988, 88, 4732-4738. (16) Roberts, A. J. SIAM J. Math. Anal. 1989, 20, 1447-1458. (17) Saaty, T. L. Modern Nonlinear Equations; Dover: New York, 1981; Chapter 5. (18) Rather than constrain the class of equations at the outset, we can always modify the equations by a judicious, if somewhat ad hoc, application of Heaviside functions, as was done in ref 2. A more elegant general solution is sought here. (19) Note that τmax may be infinite in the case of distributed delays. (20) Brown, A. J. J. Chem. Soc. 1902, 81, 373-388. (21) Henri, V. C. R. Hebd. Seances Acad. Sci. 1902, 135, 916-919. (22) Okuda, M. Prog. Theor. Phys. 1982, 68, 1827-1840. (23) Nguyen, A. H.; Fraser, S. J. J. Chem. Phys. 1989, 91, 186-193. (24) Roussel, M. R.; Fraser, S. J. J. Phys. Chem. 1991, 95, 8762-8770. (25) Roussel, M. R. Ph.D. Thesis, University of Toronto, 1994. (26) Ouellet, L.; Laidler, K. J. Can. J. Chem. 1956, 34, 146-150. (27) Darvey, I. G. J. Theor. Biol. 1968, 19, 215-231. (28) Ryde-Pettersson, U. Eur. J. Biochem. 1989, 186, 145-148. (29) Roussel, M. R.; Fraser, S. J. J. Chem. Phys. 1990, 93, 1072-1081. (30) Roussel, M. R.; Fraser, S. J. J. Chem. Phys. 1991, 94, 7106-7113. (31) A very simple strategy was adopted for numerical integration: The integrator was a straightforward constant-step-size implicit trapezoidal code. Values obtained by integration were stored in a circular buffer. The routine which returns the function to be integrated extracts the history-dependent terms from the buffer by linear interpolation. To verify the plausibility of the trajectories obtained, the results of the integration were compared to analytic results whenever possible. (32) Fraser, S. J.; Roussel, M. R. Can. J. Chem. 1994, 72, 800-812. (33) These are functions with a random component rather than truly random functions. These functions were constructed in the following manner: Since a DDE integrator uses stored points and interpolation to compute the history-dependent terms, it is necessary only to store a set of points to initialize the integrator. Random initial values for the concentrations are chosen in some positive interval. These values are then incremented by random amounts uniformly selected in a small interval at each mesh point. Negative values, when they arise, are replaced by small positive values in order to maintain the nonnegativity constraint imposed on these models. Finally, the stored values are scaled to satisfy any additional normalization conditions. (34) This procedure is similar to numerical methods for constructing the slow manifold of enzymic mechanisms; see Roussel, M. R. M.Sc. Thesis, University of Toronto, 1990. It is also closely related to seminumerical methods for constructing integral manifolds of PDE systems. The boundary of the manifold is constructed (analogous to finding an initial function which generates the slow manifold), and the manifold is then found by integration from the boundary; see Constantin, P.; Foias, C.; Nicolaenko, B.; Temam, R. Integral Manifolds and Inertial Manifolds for DissipatiVe Partial Differential Equations; Springer-Verlag: New York, 1989. (35) Wischert, W.; Wunderlin, A.; Pelster, A.; Olivier, M.; Groslambert, J. Phys. ReV. E 1994, 49, 203-219. (36) Field, R. J.; Noyes, R. M. J. Chem. Phys. 1974, 60, 1877-1884. (37) Field, R. J.; Fo¨rsterling, H.-D. J. Phys. Chem. 1986, 90, 54005407. (38) Luo, Y.; Epstein, I. R. AdV. Chem. Phys. 1990, 79, 269-299.
8330 J. Phys. Chem., Vol. 100, No. 20, 1996 (39) This scaling was chosen to make the steady-state values of the concentration variables of unit magnitude. (40) Because of the extreme stiffness of the Oregonator and delayed Oregonator equations, a slightly different numerical integration strategy had to be adopted than for the other models studied here. A third-order Gear method was used with a conservative step-size adjustment policy: The step size was decreased whenever either the solution threatened to become negative or the functional iteration routine at the core of the integrator failed to converge reasonably quickly. After a suitably large number of steps, the step size was increased again. Otherwise, the method was exactly as described in note 31 above. (41) Busenberg, S. N.; Travis, C. C. J. Math. Anal. Appl. 1982, 89, 4666. (42) Cooke, K. L.; Grossman, Z. J. Math. Anal. Appl. 1982, 86, 592627.
Roussel (43) Macdonald, N. Biological Delay Systems: Linear Stability Theory; Cambridge: Cambridge, 1989; Chapter 10. (44) The universal acceptance of this idea is relatively modern. In the early days of enzymology, it was believed by many that enzymes catalyzed only exothermic reactions. For a discussion, see Fruton, J. S. Molecules and Life; Wiley: New York, 1972; p 83. (45) Noyes, R. M. J. Phys. Chem. 1977, 81, 2315-2316. (46) Hearon, J. Z. Ann. N.Y. Acad. Sci. 1963, 108, 36-68. (47) Shear, D. J. Theor. Biol. 1967, 16, 212-228. Higgins, J. J. Theor. Biol. 1968, 21, 293-304. (48) Shear, D. B. J. Chem. Phys. 1968, 48, 4144-4147. (49) Be´lair, J. J. Dyn. Diff. Eq. 1993, 5, 607-623.
JP9600672
