Perturbative-Iterative Computation of Inertial Manifolds of Systems of Delay-Differential Equations with Small Delays

: Delay-differential equations belong to the class of inﬁnite-dimensional dynamical systems. However, it is often observed that the solutions are rapidly attracted to smooth manifolds embedded in the ﬁnite-dimensional state space, called inertial manifolds. The computation of an inertial manifold yields an ordinary differential equation (ODE) model representing the long-term dynamics of the system. Note in particular that any attractors must be embedded in the inertial manifold when one exists, therefore reducing the study of these attractors to the ODE context, for which methods of analysis are well developed. This contribution presents a study of a previously developed method for constructing inertial manifolds based on an expansion of the delayed term in small powers of the delay, and subsequent solution of the invariance equation by the Fraser functional iteration method. The combined perturbative-iterative method is applied to several variations of a model for the expression of an inducible enzyme, where the delay represents the time required to transcribe messenger RNA and to translate that RNA into the protein. It is shown that inertial manifolds of different dimensions can be computed. Qualitatively correct inertial manifolds are obtained. Among other things, the dynamics conﬁned to computed inertial manifolds display Andronov–Hopf bifurcations at similar parameter values as the original DDE model.


Introduction
Differential-equation models for the dynamics of chemical or biochemical systems often display a hierarchy of relaxation times [1][2][3]. Geometrically, this hierarchy of relaxation times corresponds to a progressive dimensional collapse, in which the solutions are confined to manifolds in phase space of lower and lower dimensions [1,2]. If we are only interested in processes occurring on longer time scales, and not in the initial fast transients, then a model constrained to a low-dimensional manifold allows us to study the important dynamics in a reduced model. Additionally, these slow-timescale models are less stiff than their parent models [4].
The "slow manifolds" [5] referred to above are invariant manifolds of the differential equations. An invariant manifold is a smooth geometric object that is mapped into itself by the action of the evolution equations [6] (Chapter 6), for instance a set of ordinary differential equations (ODEs) modeling the kinetics of a well-mixed chemical system. The simplest case of an invariant manifold is a complete trajectory integrated forward and backward in time from some initial condition. A two-dimensional invariant manifold would be generated by integrating forward and backward in time from a smooth curve transversal to the flow [6] (Chapter 6). For illustrative purposes, let us consider a three-dimensional system. Suppose that, for some set of parameters, there is a globally stable equilibrium point with one fast and two slow relaxation modes. Near the equilibrium point, these modes can be determined from the linearization of the ODEs. The fast mode will drive the system toward a two-dimensional slow invariant manifold along which the equilibrium point will be approached. If we can determine the coordinates of this manifold, then the dynamics can be reduced to a two-dimensional model capturing the slow time evolution of the original three-dimensional model after decay of the fastest transient [7].
The two-dimensional manifold may, in turn, contain a slow and a slower mode. If so, we would observe the collapse of the trajectories onto a one-dimensional manifold associated with the slower mode, and we could obtain a one-dimensional reduced model by computing the coordinates of the manifold, and then develop a time-evolution equation for motion along the manifold.
In many systems of ODEs, the stability of the equilibrium point depends on some parameters. Loss of stability when a parameter crosses a bifurcation point is associated with a finite, and usually small, number of modes. Let us say that an n-dimensional system has d unstable modes. The n − d stable modes will then drive the system to a d-dimensional unstable manifold, i.e., an invariant manifold containing the unstable modes. Since motion along the stable directions is annihilated as trajectories collapse onto the unstable manifold, motion in this manifold governs the long-term behavior of the system and we can again think about obtaining a d-dimensional reduced model for the motion confined to the manifold.
For example, in an Andronov-Hopf bifurcation, a stable focus loses stability and a stable limit cycle is born [8]. On one side of the bifurcation, the oscillatory transient approach to the stable equilibrium point occurs in a two-dimensional slow manifold; at the bifurcation itself, where the equilibrium point loses stability, the nontrivial dynamics is again in the two-dimensional manifold, which is now called a center manifold [9]; and on the other side of the bifurcation, the limit cycle is contained in the same two-dimensional manifold, now called an unstable manifold. The model reduced to the two-dimensional manifold thus captures the bifurcation behavior of the model. The manifold itself depends on the model parameters and may deform as we scan those parameters, but barring a bifurcation that changes the hierarchy of time scales of the model, the manifold persists and continues to form a basis for a reduced model [10][11][12].
In gene expression modeling, it is often necessary to include delays to take into account the nontrivial amount of time required to transcribe and translate a gene [13][14][15][16][17][18]. Delay-differential equations (DDEs), which are often used to model such systems, are instances of infinite-dimensional dynamical systems, best conceptualized as inducing a time evolution in a function space [6,19]. A DDE model would nevertheless be formulated in terms of a finite number of variables defining a state space, typically the concentrations of various chemical species in a gene expression model. It is often observed that the behavior of a DDE model, after decay of transients, lies in a low-dimensional manifold in the state space. In some situations, such as a system with an Andronov-Hopf bifurcation, this is in fact guaranteed, at least in the neighborhood of the bifurcation [20]. A differentiable state-space manifold to which the long-term behavior of DDE solutions are confined is called an inertial manifold [21]. The condition for the existence of an inertial manifold in an infinite-dimensional dynamical system is essentially the same as the condition for the existence of a slow manifold in a finite-dimensional system [22]: there must be a separation of time scales between fast relaxation and the modes governing the long-term evolution [23].
It will be briefly noted that the literature is not entirely consistent with respect to the use of the term "inertial manifold". Many authors use "invariant manifold" when describing inertial manifolds of DDEs. For state-space manifolds of DDEs, "invariance" is a problematic concept because the time-evolution operator acts on a function space [6,19]. The latter is therefore the natural setting for a concept of invariance. "Inertial manifold" is thus strongly to be preferred when describing manifolds in the state space of a DDE. On a related note, when there is no opportunity for confusion, the adjectives "invariant" or "inertial" may be dropped. Thus, in the rest of the paper, "slow manifold" may be read as shorthand for either "slow invariant manifold" or "slow inertial manifold", as appropriate, and similarly for center and unstable manifolds.
Constructing an inertial manifold for a system of DDEs is not a trivial undertaking, largely because the phase space, which is a function space, and state space, typically R n , are very different in nature. Much of the past work in this area has focused on center manifolds, where the existence of purely imaginary eigenvalues of the linearization can be exploited to bootstrap a series solution [24,25]. Qesmi and coworkers [26] and Campbell [27] have described the automation of the relevant tasks using the symbolic algebra system Maple. Farkas has proven that, in the case of a DDE with an unstable manifold, numerical integration eventually reaches this manifold [28], which opens up methods for constructing unstable manifolds based on integration, for instance the numerical projection method of Gear and coworkers [29]. Krauskopf and Green developed a method for computing one-dimensional unstable manifolds of delay-differential equations [30], which was subsequently added to the DDE bifurcation analysis package DDE-BIFTOOL [31]. Sahai and Vladimirsky have described methods based on integration to compute the two-dimensional unstable inertial manifolds of DDEs [32]. Chicone has constructed inertial manifolds as singular perturbation series, using the delay as the singular perturbation parameter [33]. These are just a handful of techniques and applications. With some exceptions, they tend to be highly specialized to one type of manifold, and sometimes to manifolds of a specific dimension. This does not negate their usefulness, but it leaves open the question of whether a more general approach might be available.
The problem of computing inertial manifolds of DDEs would be considerably facilitated if an approximating ODE could be constructed. This would open up the very large bag of tricks developed for ODE systems. One very interesting recent line of research in this direction is the Galerkin approximation method of Chekroun and coworkers [34], but this is not the path we will be following today. Rather, we will return to an old idea [35]: Assuming that the delay is not too large, we can take a Taylor expansion of a DDE in powers of the delay. This immediately reduces the DDE to an ODE. This procedure is not without its risks. Indeed, if we directly integrate such an equation for any order of expansion beyond first order in the delay, we will eventually get spurious "runaway" solutions [36][37][38]. Since truncation of the Taylor series after the first-order term was known not to generate spurious solutions, this lowest-order approximation of the DDEs has been applied to compute inertial manifolds [32,39]. However, El'sgol'ts' famous result [36] only applies to directly integrating the approximating ODEs. Chicone and coworkers have shown that expanding the delayed term in a series in powers of the delay results in correct approximations of the state-space vector field implied by a DDE [33,38,40]. Accordingly, methods for the construction of invariant manifolds that do not rely on integration should generally work with higher-order expansions.
Many gene expression models are essentially chemical models with delayed terms. Many excellent reviews of methods for computing invariant manifolds tailored to chemical models have been published [41][42][43][44][45][46][47]. The method used here and in a previous study [39] is based on Fraser's functional iteration method [1,5]. In favorable cases, this method tends to converge very quickly. Slow convergence can often be accelerated [48], and divergence can also be addressed using straightforward stabilization methods [49,50]. The approach to constructing inertial manifolds taken here can thus be thought of as combining a perturbative treatment of the delay with an iterative method for computing invariant manifolds.
As the foregoing text suggests, computing inertial manifolds of DDEs is a technically challenging undertaking. The widespread availability of DDE integrators and of sophisticated analysis tools such as DDE-BIFTOOL [51] might be thought to make an inertial manifold analysis redundant. However, there is a great deal of insight to be gained from knowing the shape of an inertial manifold, particularly when the manifold is of sufficiently low dimension to be easily visualized. In some cases, readily available experimental methods may not be able to access fast transients. Motion on the inertial manifold then becomes the experimentally observable behavior, and the geometry of the manifold may play a role in determining the stability of a steady state, for example [49]. On occasion, understanding the geometry of the slow manifold in a particular model may lead to the discovery of a class of systems expected to display similar behavior because of similar slow manifold geometries [52]. In other cases, understanding the geometry of slow manifolds may support the design of experiments to distinguish between mechanisms [53]. Additional arguments for the utility of invariant manifold calculations have been presented by Roberts [7].
In order to set the stage for the application to DDEs, Section 2 reviews Fraser's iterative method for approximating slow manifolds. The model that will be studied in that section is the zero-delay limit of the DDE model of an inducible enzyme used as our central example. In Section 3, two variations of the DDE model are studied, a three-variable model used to establish some basic techniques (Section 3.1), and a more realistic four-variable model (Sections 3.2 and 3.3). In the latter case, both three-and two-dimensional inertial manifolds are constructed. Whenever possible, insights into model behavior gained from the calculations will be pointed out.

Iterative Computation of Attracting Invariant Manifolds of ODEs
Before we proceed to the delay-differential equation problem, it is useful to review the iterative method of solution for attracting invariant manifolds of ODEs. A simple model for an inducible enzyme will be used as an example. Suppose that the cell produces an enzyme, E, the purpose of which is to metabolize a substrate S, but E is only made when S is present. Examples of such enzymes include β-galactosidase, which is required to metabolize lactose and is only produced when lactose is available [54,55], and Hmp, an enzyme involved in NO detoxification [56], which is only produced when NO starts to accumulate in a cell [57]. A simple model might be the following: The messenger RNA (R) that codes for E is synthesized at a rate that depends on the concentration of the substrate. The RNA has some natural decay rate, which may include the effect of dilution by growth.
In this and subsequent reactions, standard mass-action kinetics will be assumed, except for the empirical rate of transcription, v R (S). A k with an identifying subscript over an arrow denotes a mass-action rate constant. Concentrations of chemical species will be denoted by the corresponding italic letter, for example, S is the concentration of the substrate S.
The RNA is translated to the enzyme E that will be responsible for metabolizing S.
We consider a constant source of substrate, either as a byproduct of metabolism, or entering the cell from the surrounding medium.
The enzyme converts S to a metabolite P by the usual Michaelis-Menten mechanism: The enzyme is also subject to degradation and dilution by growth. We assume here that both the enzyme E and enzyme-substrate complex C have identical decay kinetics: The corresponding mass-action rate equations are The model is completed by specifying the dependence of the RNA synthesis rate on the substrate concentration. We assume Hill-function kinetics, which can arise in a variety of biochemical scenarios [58][59][60], and which commonly appears in gene expression models: Similar models have a long history in theoretical biology [61,62].
Although not strictly necessary, it is convenient to consider the total amount of enzyme instead of the free enzyme concentration, for reasons that will shortly become clear. In terms of this variable, the governing equations become We now transform these equations to dimensionless form following typical principles of scaling [6,[63][64][65] to obtain equations on the system's slow time scale. From the point of view of computing invariant manifolds, this is not strictly necessary, but it will facilitate some later arguments about appropriate representations of manifolds.
The maximum values of R and E T , assuming that the system is not "spiked" with one of these species from an external source, are as follows, obtained from the rate equations (5), and considering that the maximum value of v R (S) is V R : Thus, we define O(1) dimensionless variables Note that Equation (4) implies that c ≤ e.
There is not a single natural scale for S. For some parameters at low initial concentrations of E, S might rise to very large values before induction of enzyme synthesis returned its value to the vicinity of the steady state. The concentration of S at the steady state satisfies an equation that could be dominated by different terms, depending on the parameter values. Since K R controls the concentration at which RNA synthesis is half-maximal, i.e., the middle of the switching range from a low synthesis rate to a high synthesis rate, it defines a key concentration scale for the dynamics of this system. Accordingly, the scaling is at least a sensible choice. If we substitute these new variables into the rate equations, we get Proteins typically have much longer lifetimes than their mRNAs [66,67]. Moreover, both enzyme catalysis and transport are usually fast processes, so the time scales that could be extracted from Equation (8c) would be fast time scales. It follows that the slowest time scale will, in most cases, be k −1 dE . Accordingly, define the dimensionless time Also define the dimensionless parameters where K E = k −1 /k 1 and K S = (k −1 + k −2 )/k 1 are the classical Michaelis constants, the first in the equilibrium approximation, and the second in the steady-state approximation. The final set of dimensionless equations is where and dots represent differentiation with respect to θ. Typical time scales for these processes would be such that ε 1 ε 2 1. The magnitudes of α, β and κ on the other hand would be highly problem-specific. In particular, if β 1, the fastest time scale would be associated with βε 1 , i.e., c would be the distinguished "fast" variable. This is arguably the most common case in practice.
An aside on "fast" variables: Singular perturbation theory has proven to be a rich framework for discussing the invariant and inertial manifolds of differential equations [2,12,[68][69][70][71][72][73]. The standard setting is a set of differential equations in which some derivatives are multiplied by small parameters, as in Equation (11) above. Informally, we say that the corresponding variables are "fast" [74]. This is a bit sloppy, for a number of reasons, not the least of which are, first, that the differential equation for a chemical species can contain both fast and slow components, and second, that the fast and slow components of a reactive system can change over time as the concentrations change. However, it is a useful shorthand for "the concentration whose rate of change is multiplied by a small constant", so we retain this terminology.
Note that, from their definitions, we must have γ > α. Moreover, k in must be smaller than the maximum clearance capacity of the system, k −2 E max , or else the concentration of substrate grows without bound. Transforming the latter relationship to the dimensionless parameters, we have Tikhonov's theorem [6,75,76], one of the foundational theorems of singular perturbation theory, provides a key tool for determining the slow-manifold structure of this problem. If ε 1 ε 2 , ε 1 1 and β = O(1), system (11) is of the formẋ with x = (r, e) and z = (s, c). Tikhonov's theorem then states that, for sufficiently small ε 1 and after decay of fast transients, the solution of system (11) tends to the solution of the differential-algebraic system consisting of Equation (13a) with the algebraic equation g(x, z, 0) = 0. In other words, the solutions approach a slow manifold which could be expressed in the form z = Z (x) [for this problem: (s, c) = (S(r, e), C(r, e))] such that Z (x) solves g(x, Z (x), 0) = 0. (We assume uniqueness of the solution, failing which nontrivial behavior may arise, such as canard orbits [77,78].) One of Fenichel's theorems shows that there exists a neighborhood of ε 1 = 0 within which this manifold persists, with a graph that depends smoothly on ε 1 [11]. Note that the decision to use E T as a dynamical variable rather than E had the result of cleanly separating fast and slow variables, the variables s and c being involved in fast chemical processes, while the variables r and e are involved in slower synthesis processes. Following this train of thought, for small ε 1 , the fast variables should lie in a manifold computed by solvingṡ =ċ = 0. Solving these equations for the fast variables c and s, we find the following: Technically, Tikhonov's theorem requires that we set ε 1 = 0 also in Equations (14) [12,79]. However, if we do so, we obtain a manifold that does not pass through the equilibrium point for finite values of ε 1 . Moreover, as ε 1 → 0, the singular manifold (14) converges to the correct limit. For our purposes, having an approximate manifold that passes through the equilibrium point at finite values of ε 1 is a key property, so we retain these terms. Figure 1 shows several trajectories of the system (11) along with s 0 (r, e). The parameters are not those of any particular system, and were chosen mostly to provide a clear visualization of the dynamics. It seems clear from panel (a) that the trajectories are attracted to a surface near s 0 , the singular limit of the two-dimensional slow manifold. Visually, it appears that s 0 is a good approximation of the slow manifold. However, if we examine panel (b), which shows the trajectory and s 0 in an e × s projection, we see that the trajectory repeatedly crosses s 0 . The reason for this is that the correct slow manifold depends on r as well as e, unlike s 0 . Similarly, the slow manifold does not correspond to a fixed value of c, as implied by Equation (14a). Given that the singular limit is at least approximately correct, it provides some information about the expected geometry of the slow manifold. Note in particular that s 0 (r, e) has a pole at e = κ(γ − α + βε 1 ) −1 . Fenichel's theorem guarantees that there exists a neighborhood of ε 1 = 0 in which the slow manifold continues to exist, and in which the graph of the slow manifold depends smoothly on ε 1 [11]. Moreover, if we develop a geometric singular perturbation series for the slow manifold [12], the leading term would be (s (0) , c (0) ) = lim ε 1 →0 (s 0 , c 0 ), with s (0) having a pole at e = κ(γ − α) −1 . Thus we expect the slow manifold for ε 1 > 0 to have a singularity near e = κ(γ − α) −1 .
(The arguments of the ratesṡ,ṙ,ė andċ are given here to emphasize the structure of the problem, and will henceforth be omitted.) This invariance equation will be satisfied by every two-dimensional invariant manifold, and in particular by the slow invariant manifold. The time derivatives in the invariance equation are given by the evolution equations (11). Conceptually, let us think of Equation (15a) as a functional equation that can be solved for S(r, e) if we know C(r, e) as well as ∂S /∂r and ∂S /∂e: We can of course make an identical statement about Equation (15b): It implies a functional equation which we can, in this case, write down explicitly: (The functional Equation (16a) is, by contrast, an implicit equation because the nonlinear term ρ(s) inṙ makes it in general impossible to isolate S explicitly from Equation (15a).) Fraser's insight was that such functional equations can be solved by iteration [1,5,80]. For instance, we could iterate as follows: Chained iteration, in which we update the iterates sequentially as shown above, tends to converge faster than parallel iteration, where S k+1 would be computed from C k . It should be noted at this point that other iterative methods for the computation of attracting invariant manifolds have been developed, some of which have superior convergence properties [2,[83][84][85]. The method used here has the twin advantages of simplicity and of generating explicit approximations in favorable cases, including some to be studied later in this paper.
As an initial guess for the iterative process, we can use the singular limit (s 0 , c 0 ). We compute the necessary derivatives using central differences, except at the edges of the computational domain where greater stability is often achieved by extrapolating the derivatives from the interior [86]. For this particular problem, the sequence of iterates converges relatively quickly in most of the physically relevant parts of the phase space. Figure 2 shows the fourth iterate computed as described above. Large derivatives are known to be potential causes of divergence in this method [48,49], so the pole in the initial guess (and in the underlying manifold) creates anomalies that tend to spread into the computational domain. Accordingly, readers will note the imperfections in the computed surfaces at smaller values of e. Since these imperfections worsen and spread into the computational domain as more iterates are computed, it is preferable in these cases to stop iteration as soon as reasonable precision is obtained in the region of interest of the manifold. The shapes of these manifolds have both expected and unexpected features. For example, as e increases, s decreases because high concentrations of enzyme efficiently remove the substrate. Also ∂s/∂r > 0 because of the direct coupling of r to s through the control of transcription by S, reaction (1a). Neither of these properties is particularly surprising. It is however much less obvious that c should be a decreasing function of e at large r, and a non-monotone function of e at small r.
Recalling that e is a scaled version of E T , the total amount of enzyme, one might expect c to be an increasing function of e: more enzyme, more enzyme-substrate complex. However, s rises sharply at small e, resulting in a much larger fraction of the enzyme being bound to substrate, and thus explaining why, at smaller values of e and small r, c increases with decreasing e. Similarly, the increase of s with r explains why c increases with r at fixed e. For moderate values of r, this increase and the concomitant increase in s are sufficient to make c, unintuitively, a decreasing function of e.
Once we have computed the manifold, we have a reduced model since, from Equation (11c,d), Figure 3 compares trajectories obtained for the exact and reduced models. Given the small number of iterates and the difficulties in computing the manifold at its steepest part, the agreement between the two trajectories is excellent.
Finally, note that it should be possible to avoid the difficulties caused by the pole in the slow manifold by transforming to a different coordinate system [86], perhaps one with a component along s 0 (e), and one locally perpendicular to this curve.

Computation of Inertial Manifolds for DDEs
The model of the last section is not entirely realistic as it assumes that RNA and protein molecules are synthesized instantaneously [reactions (1a,b)]. From the time that a substrate is sensed by the cell's molecular machinery to the synthesis of an RNA molecule, there is a delay associated in part with transcription of the gene. Transcription can take anywhere from several tens of seconds in prokaryotes, up to hours for some of the larger genes in eukaryotes. Similarly, translation of a messenger RNA into a protein takes several tens of seconds. Thus, a better model than (5) would be with transcription and translation delays τ 1 and τ 2 , respectively. Using the same transformations and parameter definitions as before, but replacing (7a) by and introducing the total gene expression delay τ = τ 1 + τ 2 , we get the delay-differential equations where δ = k dE τ. Note that the transformation (20) has reduced the original two-delay system to a system with a single dimensionless delay, δ.

The "Fast RNA" Case
Before analyzing system (21), we consider a special limit from which we will learn some valuable lessons: Suppose that the RNA dynamics is fast, i.e., that ε 2 is vanishingly small. Then, using a standard result of singular perturbation theory for delay-differential equations [87,88], we can write the simplified system This system is almost identical to one studied previously [39], except that the older model assumed that only the free enzyme decayed. This allowed the enzyme-substrate complex (C) to act as a protected reservoir of the enzyme, so that there was no upper limit on the allowed values of k in . This is not the case for the family of models studied in this contribution, which are subject to the constraint (12).
Like the previously studied model [39], the system (22) displays an Andronov-Hopf bifurcation. Thus, it should have a two-dimensional (slow/center/unstable) inertial manifold in the vicinity of this bifurcation. Let us suppose, for the sake of argument, that βε 1 min(ε 1 , 1). Then we would expect to have an inertial manifold parameterizable in the form c = C(e, s). Our objective is to obtain an ODE model for the motion on the inertial manifold. Thus, we must eliminate Equation (22c) which contains the delayed term. Suppose that we have obtained the equation of the manifold, c = C(e, s). We can then obtain an ODE system representing the motion on this manifold by using Equations (22a,b), and solving the equation c = C(e, s) for e to close the system.
Our first task then is to obtain an approximation for the equation of the inertial manifold. We start, as usual, from the invariance equationċ In this system,ė depends on delayed values of s. This is a problem because Equation (23) is, in principle, a PDE for a manifold in the state space. Consider however expanding the nonlinearity in small powers of the delay: As previously mentioned, we cannot substitute this expansion beyond the O(δ) term into Equation (22c) and numerically integrate the equations, as shown by El'sgol'ts [36] (pp. 284-285). For this reason, previous studies had truncated the small-delay expansion after the O(δ) term [32,49]. However, we can use this expansion to approximate the vector field implied by the delay equations [38], and in particular to compute inertial manifolds of the DDEs [33,40].
The quantitiess and ... s (both evaluated at δ = 0) appearing in Equation (24) can be obtained by repeated differentiation of Equation (22a). For this particular model, we geẗ and ...
The approximate manifolds C (m) k will now have two indices: k to label the iterate, and m to indicate the last power of δ included in the expansion (24). The expansion contains powers of c, so in general it is not convenient (if it is even possible) to write down explicitly the functional equation. However, for m = 1, we have, after labeling for iteration, The m = 0 case is recovered by setting δ = 0 in this equation. If we start iteration from the initial function C (1) 0 = 0, we get the following as the first iterate: A zero initial function is often a good choice for initializing this style of iteration. In the case of ordinary chemical models, this naïve initial function generally returns the steady-state approximation as the first iterate [1,5,80], which is a useful property of this class of iterative schemes.
Note that Equation (28) is readily invertible with respect to e except when s = 0. There is therefore some hope that more accurate expressions for the inertial manifold generated by iteration will also be locally invertible, and therefore that the model reduction scheme described above will work.
Equation (27) can be iterated analytically using a symbolic algebra language. The m = 2 functional equation gives a quadratic equation in c. It could also be iterated analytically, but numerical iteration was used for m = 2 and m = 3 instead. Iteration of Equation (27) and of its implicit counterparts for m > 1 converges very quickly. The differences between the manifolds obtained for different values of m are essentially negligible (Figure 4).
We can draw a few lessons from these calculations: 1. Because the manifolds computed for different values of m are essentially identical, it appears that the higher-order corrections to the vector field are unnecessary to capture the dynamics, at least in the parameter regime of Figure 4.

2.
The calculation for m = 1 (even if carried out numerically) is much better behaved than for larger values of m. In the numerical calculations with larger m, localized imperfections appear in the manifolds. A simple routine was developed for removing these imperfections, essentially by interpolating through the region where they appear, but we can see from Figure 4b that there are small residual artifacts in the surfaces. (C 3 , on the other hand, is perfectly smooth.) 3.
If we integrate the model reduced to the inertial manifold, the m = 1 reduced model best reproduces the trajectories of the three-variable DDE model ( Figure 5). After an unphysical excursion to negative concentrations, the m = 1 model reproduces the trajectory of the three-variable model reasonably well. On the other hand, the m = 2 model has a limit cycle, which does not appear in the three-variable model until δ is above 0.54, and the m = 3 model is excessively damped. The negative concentrations in the reduced model may be due to a poor representation of the inertial manifold near the c = 0 axis, and would be an interesting subject for future investigations. The more serious misbehavior of the m = 2 and m = 3 models is likely due to a compounding of numerical issues, mostly the difficulties in solving the functional equations accurately in certain parts of the surface, numerical error during the inversion of the relation k (e, s), and of course the interaction of these two factors.
Finally, note that the δ dependence of the inertial manifold is dynamically significant. The model without a delay (i.e., m = 0) does not have an Andronov-Hopf bifurcation, but the m = 1 ODE model does, as seen in Figure 6. The limit cycle in the m = 1 model is a little larger than that in the DDE model, and is a little distorted. Nevertheless, there is clearly a wide range of parameters near the Andronov-Hopf bifurcation where the two models have similar behavior. In past work, it was found that an m = 1 model accurately reproduced the Andronov-Hopf bifurcation curve of a related DDE model [39]. Although such a calculation has not been done for this model in order to focus on other issues, it seems clear that this will be true for this model as well.   Figure 4. All trajectories were started from initial conditions near the origin, s(0) = 0.02, c(0) = 0. For the three-variable model, e(0) = 0 and the initial function was s(θ) = 0 for θ < 0 (i.e. there is a small discontinuity in the initial function at θ = 0). Figure 6. Trajectories of (a) the three-variable DDE model and of (b) the m = 0 and m = 1 reduced ODE models. The m = 1 model has an Andronov-Hopf bifurcation and can generate limit-cycle oscillations, like the three-variable model, whereas the m = 0 model never undergoes an Andronov-Hopf bifurcation. All parameters are as in Figure 4, except δ = 0.7.

Elimination of the Fastest Relaxation Mode
We now return to the original four-variable model including RNA, enzyme, substrate and enzyme-substrate complex, Equation (21). Let us suppose that we want to develop a state-space model, i.e., carry out a minimal reduction in which we eliminate the delay and obtain ODEs for the dynamics in a three-dimensional state space by computing a three-dimensional inertial manifold in the four-dimensional state space. As discussed previously, for "normal" parameters, i.e., parameter sets with ε 1 ε 2 1, the system has three distinct time scales: fast enzyme kinetics, moderately fast mRNA dynamics, and slow dynamics of the total amount of enzyme. In this parameter range, the fastest variable, which is the one whose geometry will be most appropriate to this reduction, will be either s or c, depending on the values of β and γ. (Either β < 1 or γ α would tend to make c the fastest variable.) However, the delay appears in the equation for r, so this is the variable in which we need to approximate the dynamics on the slow manifold in order to obtain the desired reduction. Let us say, for the sake of argument, that c is the fastest variable. We can compute an (approximate) inertial manifold in the form c = C(e, r, s) as we did in the previous section. We could then, in principle, have a reduced model that consists of the following set of differential-algebraic equations: As we did for the three-variable model, the last of these equations could in principle be numerically inverted to give us a value of r for any given c, e and s, thus closing the system. The problem is that c depends only indirectly on r through e. This has the result of making ∂c/∂r on the manifold small and difficult to compute accurately. The inversion process is therefore numerically ill-conditioned, so that a useful reduction is not obtained. Accordingly, in the usual parameter regime, a three-variable ODE model is not available.
The assumption that ε 2 is larger than ε 1 will be true for the vast majority of enzyme synthesis/catalysis systems, but is not necessarily so. Values of ε 2 in a mouse cell line can be calculated directly from the data of Schwanhausser and coworkers [67], yielding values in the range of 1.5 × 10 −3 to 46, with a median of 0.20. The same data set gives a maximum value of k dE of 3.9 × 10 −4 s −1 . For prokaryotes, values larger than this by an order of magnitude are possible [89]. Rate constants for enzyme-substrate association (k 1 ) can be as low as 10 3 M −1 s −1 [90,91]. Since enzyme concentrations in excess of 1 µM are routinely seen in vivo [92], values of ε 1 > ε 2 are certainly within the realm of possibility. In the extreme, this parameter regime leads to the three-variable model. Since ε 2 is never vanishingly small, it is worth reexamining this case using the manifold construction techniques we have been studying to compute an ODE state-space model valid when ε 2 is small, but not vanishingly so. If ε 2 is smaller than both ε 1 and βε 1 , then r is the fastest variable in the system, and we should be able to directly eliminate it by computing an inertial manifold of the form r = R(e, s, c). This ansatz leads to the invariance equationṙ The general plan is as before: we will use the expansion (24) to replace the delayed term in Equation (21c) prior to substitution into the functional equation, and then iteratively solve the invariance equation for R(e, s, c) from some suitable initial guess. If we expand ρ to first order in δ, the resulting iteration is To go beyond first order, we need and ...
Up to third order in δ, it is possible to analytically isolate R from Equation (30). The functional equations for R (m) become more and more complicated, but the work is easily carried out in a symbolic algebra system. Iteration of these functional equations can thus be carried out analytically, also in a symbolic algebra system. The size of the algebraic expressions generated by the iteration process grows rapidly, limiting the number of iterates that can practically be computed. Maple, the symbolic algebra system used in this study, lacks a function that gives a useful measure of the size of an expression. A simple workaround is to have Maple generate (say) Matlab code, and to count the number of bytes in these files. The optimize option was not used in order to measure the size of the expressions as stored by Maple. Also note that it is not feasible to use simplify() on the expression generated during iteration as the computation time required to execute this function grows very rapidly with the size of the expression. Table 1 shows the growth of the size of the expressions for the iterates with both k and m. Note particularly the enormous sizes of the expressions generated for m = 2 and 3 after a modest number of iterations. This is another reason to prefer the lowest-m expansion that suffices to reproduce the behavior of the original system.
If we take as an initial guess R (1) 0 = 0, then from Equation (31) we get the following O(δ) first approximation to the inertial manifold: This is the simplest approximation to the inertial manifold that takes account of the delay. For m = 1, there is no significant change in the shape of the manifold for k > 2. This is difficult to show directly since we cannot plot a function of three variables. We could take cuts through this function at, say, different values of c, but these are not very enlightening. Instead, we compute trajectories of the reduced system using r = R (m) k (e, s, c) for different values of k, and then calculate the scaled difference between iterates, the iteration error (IE): In this equation, (e(θ), s(θ), c(θ)) (m) k refers to a trajectory of the reduced model. This is a very stringent test because differences in frequency can cause large differences in the trajectories even if the overall dynamics are very similar. Figure 7 shows these iteration errors for a parameter set in the regime where we expect to be able to carry out this reduction. Near the initial condition, it is difficult to measure IE (m) k reliably because of the small denominator in Equation (35), so values of this statistic near θ = 0 should be ignored. For k = 2, we see large-amplitude oscillations, indicating that the two manifolds are sufficiently different to cause a loss of coherence. When we calculate the next iterate, the iteration error becomes much smaller and loses any obvious periodicity, which is evidence of rapid convergence of the iterations. The same behavior is observed for higher values of m. While iteration converges rapidly, it turns out that the reduced models are only qualitatively correct. Figure 8 shows trajectories of the DDE model as well as of the state-space ODE models. For the parameters chosen, the DDE model has a stable focus. On the other hand, the reduced models for m = 1 and 2 have limit-cycle attractors, i.e., they undergo an Andronov-Hopf bifurcation at a smaller value of δ (keeping the other parameters constant) than the DDE model, near δ = 0.84 for the m = 1 reduced model vs. δ = 1.04 for the DDE model. As for the m = 3 reduced model, it is much more strongly damped than the DDE model, and it fails completely if we try to increase δ too much.
Note that the value of ε 2 chosen for the calculations of Figures 7 and 8 is not particularly small compared to ε 1 . The m = 1 and m = 2 state-space models are qualitatively correct, even if there are some quantitative discrepancies such as the exact parameter values at the Andronov-Hopf bifurcation. The agreement improves if we decrease ε 2 .

Computing a Two-Dimensional Inertial Manifold
Given the Andronov-Hopf bifurcation in this system, the lowest-dimensional description of the dynamics that captures the behavior to either side of the bifurcation is motion on a two-dimensional inertial manifold. As argued previously, the most likely parameter regime is one in which c and s are the fast variables. However, as in the three-variable model, this leads to an implicit manifold representation problem which will be left for future study. Rather, we consider the regime where β 1 and ε 2 {ε 1 , 1}. This is the easiest biochemically plausible case since r is then one of the fast variables, which means that the delay can be directly eliminated. The other fast variable in this case is of course c. These considerations imply a parameterization of the manifold c = C(e, s), r = R(e, s). The corresponding invariance equations arė The delayed term does not appear in Equation (36a), so the functional iteration for C(e, s) is the same regardless of the delay expansion order: (C (m) k will eventually depend on m because of its dependence on R (m) k .) The functional equations for R(e, s) can also be derived analytically, although for m > 1, the derivation is best left to a symbolic algebra system. For m = 1, we have assuming chained iteration, C first. Even starting with our usual naïve initial function of R (m) 0 = 0, we see that R (1) 1 will depend on δ. Corrections to the slow manifold due to the delay therefore enter at a level equivalent to the application of the steady-state approximation.
From the functional equations (37) and (38), we can see that at this level of approximation of the vector field, corrections to the manifold due to the delay depend on δ/ε 1 . Thus, the delay will become less significant to the slow dynamics at small values of δ or large values of ε 1 . Here we see an example where the analysis tells us something that might have been difficult to intuit: In terms of the original dimensional parameters, the condition δ/ε 1 1 becomes k 1 E max τ −1 . Thus, if the association of the enzyme and substrate is sufficiently slow compared to the gene expression cycle, the delay associated with this cycle should have a negligible effect on the geometry of the slow manifold and, at least from the perspective of the motion on the two-dimensional slow manifold, an ODE model neglecting the transcription and translation delays would be appropriate. Note that this says nothing about the transients, which could still be very different between the models with and without delays.
The reduced model obtained from Equation (21) for a computed manifold (C(e, s), R(e, s)) is the planar system Given explicit functional equations, we could try to carry out the iterations analytically in a symbolic algebra system, but again the expressions quickly grow to the point where they cannot practically be handled, so numerical iteration was used. The second iterate was obtained analytically for each m and used to verify reasonable accuracy of the numerical scheme. These differences (not shown) are generally very small (O(10 −6 ) or less), except at the s = 0 edge. Note that the boundary of the computational domain is often a source of instability in this class of iterative methods [50,86]. Now consider the first iterate of Equation (38) starting from the naïve initial function: For n > 1, dρ/ds → 0 as s → 0. Thus, R 1 itself tends to zero as s → 0. Estimating the derivative in this region therefore involves taking differences of small quantities, which is difficult to do accurately, resulting in larger errors in the numerical iteration process in this region. Fortunately, these errors do not destabilize the iterative scheme, at least if we stop iteration as soon as reasonable convergence is achieved. Convergence is (again) very fast, so three or four iterates are sufficient to reduce the iteration error below the threshold where it has any observable effect on the dynamics. Because the approximated manifolds are now paired surfaces C(e, s) and R(e, s), they are readily visualized. Figure 9 shows the computed inertial manifolds at various levels of truncation, m. Despite its coupling to R (Equation (37)), C depends very little on m, so a representative approximation is shown. On the other hand, the shapes of the R (m) surfaces are dramatically different. This is a possible sign of a divergent series in δ. (Note that δ is not particularly small here.) However, which representation of the inertial manifold is the best? A simple way to answer the latter question is to look at trajectories. For the parameters of Figure 9, the DDE model has a limit cycle ( Figure 10). The system reduced to the m = 1 approximate inertial manifold also has a limit cycle (not shown), but this limit cycle wanders out of the positive quadrant of the (e, s) plane. We might accept unphysical excursions of trajectories started from some initial conditions, but an unphysical attractor is not acceptable in any potential application. The m = 2 attractor is, on the other hand, an entirely sensible-looking limit-cycle entirely contained in the positive quadrant, although it is too large and of a different shape than the limit cycle of the DDE model ( Figure 10). The significantly larger size of the limit cycle in the m = 2 reduced model is associated with an Andronov-Hopf bifurcation that occurs at a smaller value of δ than in the DDE model. For m = 3, all we get is a strongly damped oscillation (not shown). It thus seems that the m = 2 approximation is the best representation of the inertial manifold available in the context of this method for these parameters. Figure 10. Trajectories of the DDE model (21) and of the system reduced to the inertial manifold for m = 2. The parameters for the first two curves are those of Figure 9. In particular, δ = 3. For the third curve, the DDE model was integrated with δ increased to 7.
The observation that the Andronov-Hopf bifurcation is advanced in the m = 2 reduced model relative to the original DDE model raises an interesting question, particularly since the limit cycle in the reduced model has a different shape than that in the DDE model: are the dynamics of the reduced model simply shifted along the δ axis in the parameter space? Along with the two limit cycles already discussed, Figure 10 shows the attractor of the DDE model at a larger value of δ. This limit cycle has a shape much more like that of the reduced model. Remarkably then, the small-delay approximation captures the rough shape of the vector field on the inertial manifold, but at a somewhat different value of δ than the one used in constructing the manifold. This suggests that there may exist a transformation of the δ axis that would bring the two models into closer correspondence. It is not clear whether a systematic method could be developed to estimate the mapping from the true value of δ to the effective δ of the reduced model. This would be an interesting avenue for future investigation.
Finally, note that the graph of R(e, s) for the best representation of the slow manifold, R (2) (Figure 9c) has a complex shape that could not be predicted from intuition alone. The question of whether this surface could be systematically explored in realistic experiments is also left for future investigation.

Discussion
Throughout the examples considered in this paper, we consistently obtained rapid convergence of the iterative schemes at fixed m. On the other hand, the expansion of the delayed nonlinearity in powers of δ was more delicate. Sometimes, the best results were obtained for m = 1, and sometimes for m = 2. It seems likely, given the results of Chicone [33,38,40] that there is no fundamental limitation on the order of expansion in the delay, and that the "best" order of expansion therefore depends on the model and parameters. In Chicone's work, it was assumed that the delay was sufficiently small for the Taylor series to converge. Here, relatively large values of the delays have been considered, and convergence is not assured, but useful results can still be obtained provided we pick the "correct" order of expansion. The situation appears to be (at least) analogous to asymptotic series, where we get sensible approximations provided we cut off the expansion as soon as it starts to diverge [93]. Now that the method has been shown to give reasonable results in a number of challenging cases, the stage is set for a more deliberate analysis of its mathematical properties.
Despite the relatively large values of the delay used in this study, it was possible to calculate manifolds on which the dynamics are at least qualitatively correct. Note that the calculations were carried out to either side of the Andronov-Hopf bifurcation, even though results were not reported for every case tried. Thus, the method proposed here is able to represent inertial manifolds of DDEs over a range of parameters. This raises the question of whether there is some other development of the delay term that could be used. The Galerkin approximation method of Chekroun and coworkers [34] has already been mentioned. As a general point, it would be useful to have a variety of approximate transformations of delay equations valid in different parameter regimes.
Scaling played an important role throughout the work reported here. For every model variation analyzed, values of the singular perturbation parameters ε 1 , ε 2 and β would inform our choice of representations of the manifold. When the correct representation was chosen, rapid convergence of the iterative process was always obtained. In some calculations not reported here, choosing the "wrong" representation for the parameter range considered always led to difficulties with the iterative method, essentially because the manifolds were nearly vertical in parts of the state space. When analyzing smaller models for which scaling is practical, especially if we have some idea of the relevant parameter values, it is therefore well worth investing some time into a scaling analysis. This is, of course, not an original recommendation [63,65].
Implicitly, in selecting the representations of the manifolds in the delay model, we assumed that a version of Fenichel's geometric singular perturbation theory [12] and, in particular, of his persistence and smoothness results [11], extended to inertial manifolds. Others have made the same assumption [94]. Chicone has shown that similar persistence and smoothness results hold using the delay as the small parameter in the singular perturbation expansion [33]. In Chicone's approach, the singular limit is an ODE. The singular limit of systems such as (21) can be a set of DDEs with an algebraic constraint, or a set of ODEs with an algebraic functional constraint (e.g., ρ(s(θ − δ)) = r(θ) if ε 2 → 0). To the author's knowledge, persistence and smoothness of inertial manifolds for such systems has not been systematically studied. This is clearly an important area for future work. Given Chicone's work and other extensions of singular perturbation theory to DDEs [87,88], it seems likely that such an extension of Fenichel's results could be obtained.
Given the self-correcting nature of the iteration process, it is not necessary for the numerical iterates to be extremely accurate, provided they are smooth. Calculation of the derivatives is the main source of purely numerical instabilities. (There is also the possibility that the iteration scheme itself is unstable [48,49].) As noted earlier, the derivatives are currently computed by finite differences, with extrapolation used at the edges to avoid well-known difficulties with this class of iterative methods due to poor derivative estimates at the boundary of the computational domain [50,86]. Making the mesh tighter will make the derivative estimates more accurate in highly curved regions at the cost of potentially losing significant figures. However, the results are not highly sensitive to the computational mesh, provided we avoid making it too tight. What does matter is making sure to extend the mesh far enough to catch the full trajectory, even if this trajectory will wander into the unphysical, negative region of state space, as frequently happens to trajectories started near the boundary in these models ( Figures 5, 6 and 8). Since this was not expected at first, the computational domain was generally set up to only include the part of the positive orthant likely to be visited. When trajectories wandered out of this region, the simulations had to be stopped, and the manifold recomputed in an extended region of state space that would hopefully be sufficient to accommodate the entire trajectory. The combustion community has often faced a similar problem, i.e., a lack of a priori knowledge of the region of the very high-dimensional phase space to be visited in a simulation. Their solution has been to extend the representation of the manifold as needed in directions visited by trajectories [95]. When the inertial manifold of a DDE needs to be computed numerically, a similar method could be used, first building a representation of the manifold in a region of the state space that we expect to be visited, and dynamically growing the manifold as needed.
When the iterative process resulting from the combined perturbative-iterative method proposed here can be carried out analytically, we end up with a small ODE model that may be amenable to computer-aided analysis. For example, the Andronov-Hopf bifurcation in the reduced model (39) will occur when 1 − ∂R ∂e with the left-hand side evaluated at the steady state, provided the quadratic discriminant of the characteristic polynomial (for which there is also a simple expression) is negative. Given analytic expressions for C(e, s) and R(e, s) and their dependence on the parameters, the Andronov-Hopf bifurcation condition is an explicit condition. Since steady-state coordinates depend on the parameters, to calculate Andronov-Hopf bifurcation curves, we solve the steady-state problem associated with Equation (39) along with (41) as a system of three equations in three unknowns, namely e, s and one free parameter at fixed values of the other parameters. This is an entirely straightforward numerical process [39]. By contrast, the Andronov-Hopf bifurcation problem for a DDE involves a characteristic quasipolynomial (involving both polynomial terms and terms in e −λδ , where λ is the characteristic value) [6] (Chapter 12). The real and imaginary parts of this characteristic quasipolynomial must be separated, and the real part set to zero. Obtaining an explicit bifurcation condition takes quite a bit of work. Of course, so does generating an inertial manifold. Which one is preferable becomes a matter of how often we will be solving the Andronov-Hopf bifurcation condition, as well as whether we have other uses for the representation of the manifold. If we can only carry out the iteration numerically, it is in principle still possible to solve for Andronov-Hopf bifurcations using Equation (41), but it becomes much more laborious since the inertial manifold depends on the parameters. It will therefore be necessary to recompute the manifold (or at least a piece of it near the steady state) dynamically during the solution process. Instead of the iterative process used here, the method of invariant grids [96], which involves using a small number of points to represent a small part of the manifold, might be a more efficient way to represent and recompute the manifold in this case. An alternative might be to use continuation methods to compute the change in the manifold as parameters are changed [97]. It must be admitted though that working directly with the DDEs might be more efficient in these cases.
The model studied here was reducible to a single delay. However, gene expression models often involve multiple delays of different sizes [98][99][100][101] and will not in general be reducible to a single delay. The development of methods for calculating inertial manifolds would be particularly useful for such systems, since the complexity of linear stability analysis increases significantly with the number of delays. This is, to my knowledge, an entirely open problem area.

Conclusions
A method for computing inertial manifolds of systems of DDEs has been developed for models of the formẋ = f (x(θ), y(θ)) , (42a) where x is an (n − 1)-vector, y is a scalar, and n > d, the dimension of the manifold. Small parameters multiplying the derivatives, if they appear, can help guide the selection of manifold representations, but are not strictly necessary. This method involves a Taylor expansion of the delayed term in powers of the delay, followed by the use of an iterative method for the computation of the inertial manifold of the approximated vector field. In favorable cases, a computer algebra system can be used to obtain fully analytic approximations, with all of the attendant advantages. One notable property of the method developed here is that it can be applied to calculate inertial manifolds of any dimension. It is not limited to computing, for instance, two-dimensional center manifolds, as we saw in Section 3.2 where a three-dimensional manifold was obtained. At this time, the selection of the order of truncation of the Taylor series is based on ad hoc observations. Part of the difficulty in selecting the order of truncation is that the Taylor series is itself a complicated function of the state-space variables. (See, e.g., Equations (24)- (26)). Some cleverness will therefore be required to discover a systematic method for choosing the truncation order.
Funding: This research was funded by the Natural Sciences and Engineering Research Council of Canada.