Tikhonov-Fenichel reduction for parameterized critical manifolds with applications to chemical reaction networks
aa r X i v : . [ m a t h . D S ] M a y Tikhonov-Fenichel reduction for parameterizedcritical manifolds with applications to chemicalreaction networks
Elisenda FeliuDepartment of Mathematical Sciences, University of CopenhagenUniversitetsparken 5, 2100 Copenhagen, Denmark [email protected]
Niclas KruffLehrstuhl A f¨ur Mathematik, RWTH Aachen52056 Aachen, Germany [email protected]
Sebastian WalcherLehrstuhl A f¨ur Mathematik, RWTH Aachen52056 Aachen, Germany [email protected]
May 22, 2019
Abstract
We derive a reduction formula for singularly perturbed ordinary dif-ferential equations (in the sense of Tikhonov and Fenichel) with a knownparameterization of the critical manifold. No a priori assumptions con-cerning separation of slow and fast variables are made, or necessary.Weapply the theoretical results to chemical reaction networks with mass ac-tion kinetics admitting slow and fast reactions. For some relevant classesof such systems there exist canonical parameterizations of the variety ofstationary points, hence the theory is applicable in a natural manner. Inparticular we obtain a closed form expression for the reduced system whenthe fast subsystem admits complex balanced steady states.
MSC (2010):
Key words : Singular perturbation, critical manifold, chemical reactionnetwork, complex balancing Introduction
A fundamental result on singular perturbation reductions, due to Tikhonov andFenichel, allows to reduce the dimension of an ordinary differential equationwith a small positive parameter in the asymptotic limit when the parameter ap-proaches zero. This theorem has numerous applications in the sciences, in par-ticular to chemical and biochemical reaction networks, especially quasi-steadystate (QSS) for certain chemical species, and partial equilibrium approximation(PEA) for slow and fast reactions. However, the application of Tikhonov’s andFenichel’s theory to reaction networks may pose some computational problems,and the purpose of the present paper is to address and resolve one of theseproblems.Throughout the present work we assume that a suitable small parameter (ina system possibly depending on several parameters) has been identified, hencewe deal with a singularly perturbed ordinary differential equation. However, wedo not assume the equation to be given in separated fast and slow variables, sothe usual version of the reduction theorem is not directly applicable. In applica-tions, fast-slow variable separation is frequently not satisfied a priori and worse,there may be no explicit way to rewrite the system in fast-slow form. Generallyone may circumvent (and to some extent resolve) this problem by resorting toan “implicit” version of the reduction, which admits the critical submanifold ofphase space as an invariant set, but this approach may also encounter compu-tational feasibility problems. Given this background, we derive in the presentpaper an explicit singular perturbation reduction that is applicable whenever aparameterization of the critical manifold is known.The paper is organized as follows. After some preliminary work (mostly re-calling notions and results from the literature) we derive in Section 2 a generalformula for Tikhonov-Fenichel reduction when a (possibly local) parameteri-zation of the critical manifold is given, and also consider some special cases.We illustrate the procedure by some small examples, briefly indicating that therange of applications is not restricted to chemical reaction networks. However,this reduction formalism seems particularly useful for reaction networks whenthe partial equilibrium approximation is applicable, since in many instances va-rieties of stationary points admit a canonical parameterization. This setting isdiscussed in detail in Section 3, and our results include a closed form reductionformula for fast subsystems that are complex balanced, as well as a discussion oflinear attractivity properties of slow manifolds. To finish the paper we discusssome examples.
We consider the singular perturbation reduction of ordinary differential equa-tions, with no a priori assumption on separated (slow and fast) variables. Thuslet U ⊆ R n be open, ε >
0, and let h be a smooth function in some neigh-borhood of U × [0 , ε ), with values in R n . This defines a parameter-dependent2ystem of ordinary differential equations, viz.(1) ˙ x = h (0) ( x ) + εh (1) ( x ) + ε . . . , x ∈ U, ε ≥ , and rewritten in slow time scale τ = εt we have a singularly perturbed system(2) x ′ = 1 ε h (0) ( x ) + h (1) ( x ) + ε . . . , x ∈ U, ε ≥ . Here h (0) is called the fast part and h (1) the slow part of either system. Wefocus on the behavior of the (1), (2) as ε →
0, and we will restrict attention toscenarios for which, modulo a coordinate transformation, the classical singularperturbation theorems of Tikhonov [22] and Fenichel [9] are applicable. (To bespecific, we refer to the version of Tikhonov’s theorem as given in Verhulst [23],Theorem 8.1 ff.; see also [11], Section 2.1.) While one can establish intrinsic con-ditions for the existence of a coordinate transformation to “Tikhonov standardform” with separated slow and fast variables, such a transformation cannot gen-erally be obtained in explicit form. We first recall a general implicit reductionprocedure developed in [11, 19], and then present, as a new result, a version ofthe reduced system that can be computed explicitly when a parameterizationof the critical manifold is known.
We recall the essential results on coordinate-independent Tikhonov-Fenichel re-duction from [11, 19]; in particular we refer to [11], Theorem 1 and the subse-quent remarks.
Proposition 1.
Let system (1) be given, and denote by V ( h (0) ) the zero set of h (0) . Moreover let < r < n and set s := n − r > .(a) Assume that a ∈ V ( h (0) ) has the following properties. • There exists a neighborhood e U of a such that rank Dh (0) ( x ) = r for all x ∈ Z := V ( h (0) ) ∩ e U ; in particular Z is an s -dimensional submanifoldof R n . • For all x ∈ Z there is a direct sum decomposition R n = Ker Dh (0) ( x ) ⊕ Im Dh (0) ( x ) . • For all x ∈ Z the nonzero eigenvalues of Dh (0) ( x ) have real part < .Then in some neighborhood of a there exists an invertible coordinate trans-formation from (1) to Tikhonov standard form ˙ y = εf ( y , y ) + O ( ε )˙ y = f ( y , y ) + O ( ε ) with separated slow and fast variables; moreover the fast system satisfies alinear stability condition. b) Conversely, the conditions in part (a) are necessary for the existence of alocal coordinate transformation to Tikhonov standard form.(c) One may choose e U such that there exists a product decomposition with func-tions µ ( x ) taking values in R r × , P ( x ) taking values in R n × r , such that (3) h (0) ( x ) = P ( x ) µ ( x ) , for all x ∈ Z ; moreover rank P ( a ) = r , rank Dµ ( a ) = r and Z = V ( µ ) ∩ e U .
Here the entries of µ may be taken as any r entries of h (0) that are func-tionally independent at a .(d) The following system (in slow time) is defined on e U , and admits Z as aninvariant set: (4) x ′ = (cid:0) I n − P ( x ) A ( x ) − Dµ ( x ) (cid:1) h (1) ( x ) , with A ( x ) := Dµ ( x ) P ( x ) . The restriction of this system to Z corresponds to the reduced equation inTikhonov’s theorem. We refer to Z as the local critical manifold (or local asymptotic slow mani-fold) of system (1). Remark 1. (a) Note that Dh (0) ( x ) = P ( x ) Dµ ( x ) on Z , due to µ ( x ) = 0.Since P ( x ) has full rank on Z , Dh (0) ( x ) and P ( x ) have the same columnspace.(b) The eigenvalues of A ( x ), x ∈ Z , are the nonzero eigenvalues of Dh (0) ( x )whenever the latter has rank r ; see [11], Remark 3.(c) We call(5) Q ( x ) := I n − P ( x ) A ( x ) − Dµ ( x )the projection operator of the reduction. For each x this is a linear projectionof rank s = n − r which sends every element of R n to its kernel componentfrom the kernel-image decomposition with respect to Dh (0) ( x ).(d) Formally system (4) is defined whenever A ( x ) is invertible, and by Fenichel’sresults it corresponds to a reduced system as ε → A ( x ) have nonzero real part (normal hyperbolicity).4 emark 2. The reduced system may just have the form x ′ = 0; in particularthis occurs in the following scenario: h (0) always admits n − s independentfirst integrals near any point of Z , and locally every point of Z is uniquelydetermined as an intersection of Z with suitable level sets of these first integrals;see [11], Subsection 2.3. Now, in the special case when h (1) admits the same firstintegrals, then Z as well as every intersection of Z with level sets is invariant forthe reduced equation, meaning that every point of Z is invariant, thus stationary.However, the only information to be gained from x ′ = 0 for small ε > ε or higher. (Generally the reduced system (4) in slow time represents only the O (1) term in ε .)While Proposition 1 provides a general coordinate-free approach to singularperturbation reduction, the critical manifold Z is given only implicitly via thezeros of h (0) , and one cannot generally expect an explicit reduction to a system in R s . Moreover there may be a problem with the feasibility of the computations,in particular with the computation of the projection matrix Q . Therefore it isnatural to search for simplified reduction procedures in special circumstances.One notable scenario appears when a parameterization for the critical manifoldis explicitly known, and we will next discuss reduction in this case. We keep the assumptions and notation from Proposition 1, in particular thedecomposition (3), the s -dimensional local critical manifold Z (being the zeroset of µ , as well as of h (0) ), and the reduced system(6) x ′ = Q ( x ) h (1) ( x ) on Z. Now assume that there is an open set W ⊆ R s and a smooth parameterization(7) Φ : W → Z, rank D Φ( v ) = s for all v ∈ W. Then every solution x ( t ) of (6) with initial value in Φ( W ) can be written in theform x ( t ) = Φ( v ( t )) , for t in some neighborhood of 0, and differentiation yields(8) D Φ( v ( t )) v ′ ( t ) = x ′ ( t ) = Q (Φ( v ( t ))) · h (1) (Φ( v ( t ))) . The remaining task is to simplify this expression.
Theorem 1. (a) For every v ∈ W there exists a unique R ( v ) ∈ R s × n such that Q (Φ( v )) = D Φ( v ) · R ( v ) . (b) The reduced system, in parameterized version (8) , is given by (9) v ′ = R ( v ) · h (1) (Φ( v )) . c) The matrix R ( v ) is uniquely determined by the conditions R ( v ) · P (Φ( v )) = 0 and R ( v ) · D Φ( v ) = I s , and therefore can be obtained from the matrix equation R ( v ) · ( D Φ( v ) | P (Φ( v ))) = ( I s | with ( D Φ( v ) | P (Φ( v ))) invertible. In particular, v R ( v ) is smooth.(d) For every x ∈ Z let L ( x ) ∈ R s × n be of full rank s and such that L ( x ) Dh (0) ( x ) =0 ; equivalently L ( x ) P ( x ) = 0 . Moreover define L ∗ ( v ) := L (Φ( v )) . Then R ( v ) = (cid:0) L ∗ ( v ) D Φ( v ) (cid:1) − L ∗ ( v ) , and the reduced system, in parameterized form, is given by (10) v ′ = (cid:0) L ∗ ( v ) D Φ( v ) (cid:1) − L ∗ ( v ) h (1) (Φ( v )) . Proof.
For every v ∈ W one has h (0) (Φ( v )) = 0, and by differentiation Dh (0) (Φ( v )) D Φ( v ) = 0 . Thus the image of D Φ( v ) is contained in the kernel of Dh (0) (Φ( v )), and thesetwo vector spaces have dimension s , hence they are equal. In turn, for x ∈ Z the kernel of Dh (0) ( x ) is by construction equal to the image of Q ( x ). Thus,for every v the matrices Q (Φ( v )) and D Φ( v ) have the same column space, andthe latter has full rank. Therefore every column of Q (Φ( v )) is a unique linearcombination of the columns of D Φ( v ). Rewritten in matrix language, this is theassertion of part (a). Part (b) is now obvious from equation (8) and injectivityof D Φ( v ).The first condition given in part (c) is a consequence of part (a), the identity Q ( x ) · P ( x ) = 0 for all x ∈ Z (which is readily verified from the defining equation(5)), and the fact that D Φ( v ) is an injective linear map. The second conditionfollows from the fact that D Φ( v ) · R ( v ) = Q (Φ( v )) is a projection of rank s , byusing Lemma 1 in the Appendix. Invertibility of the matrix ( D Φ( v ) | P (Φ( v )))follows from the direct kernel–image decomposition with respect to Dh (0) (Φ( v )),since the columns of D Φ( v ) span the kernel and the columns of P (Φ( v )) spanthe image (see [11] for more details).To prove (d), first recall from Remark 1 that Dh (0) ( x ) and P ( x ) have thesame column space, therefore L ( x ) P ( x ) = 0 on Z . This and R ( v ) P (Φ( v )) = 0from part (c) imply that R ( v ) = Λ( v ) L (Φ( v )) for all v ∈ W , with Λ( v ) ∈ R s × s uniquely determined. Using now the second condition in part (c), we haveΛ( v ) L (Φ( v )) D Φ( v ) = I s , hence L (Φ( v )) D Φ( v ) is invertible andΛ( v ) = (cid:0) L (Φ( v )) D Φ( v ) (cid:1) − , which leads to the asserted expression.6 emark 3. (a) To determine the reduced equation (10), there is no need forexplicit knowledge of the projection matrix Q , or of the matrix P from thedecomposition. However, the column space of Dh (0) ( x ), x ∈ Z , is a crucialingredient.(b) On the other hand, knowledge of P and µ seems indispensable for the com-putation of A ( x ) = Dµ ( x ) P ( x ), and of A (Φ( v )). Note that the eigenvaluesof the latter provide direct information on the stability of the critical man-ifold; see Remark 1 (b).We consider some special cases in more detail. Corollary 1.
Assume that Φ( v ) = (cid:18) Φ ( v )Φ ( v ) (cid:19) , with Φ ( v ) ∈ R s and D Φ ( v ) invertible, for all v ∈ W, and partition P ( x ) = (cid:18) P ( x ) P ( x ) (cid:19) with P ( x ) ∈ R s × r . Then R ( v ) = ( R ( v ) | R ( v )) with R ( v ) = I s − P (cid:0) D Φ D Φ − P − P (cid:1) − D Φ D Φ − R ( v ) = P (cid:0) D Φ D Φ − P − P (cid:1) − . In these expressions the argument of D Φ and D Φ is v and the argument of P and P is Φ( v ) .In the special case when Φ ( v ) = v we get R ( v ) = I s − P ( D Φ P − P ) − D Φ R ( v ) = P ( D Φ P − P ) − . Proof.
With R i given as above one verifies R D Φ + R D Φ = I s and R P + R P = 0by direct computation. Rewriting, one obtains( R | R ) (cid:18) D Φ P D Φ P (cid:19) = ( I s | , and this is the defining property of R in Theorem 1. Remark 4. (a) Up to a relabeling of variables in R n , a partitioning for Φ( v )as required in Corollary 1 always exists locally, due to the rank conditionon the derivative. 7b) The special case Φ ( v ) = v occurs when the critical manifold is the graphof some function. For this case, reduction formulas were derived earlier byFenichel [9], Lemma 5.4, and Stiefenhofer [21], Equation (2.13).(c) In the yet more special case that Φ ( v ) = v and Φ ( v ) = 0 the procedureyields the familiar quasi-steady state reduction. Indeed, in this case one has µ ( x ) = x for x = ( x , x ) tr and x ∈ R s , thus R = I s , R = − P P − andthe reduced equation is v ′ = (cid:0) I s | − P (( v, tr ) P − (( v, tr ) (cid:1) h (1)1 (( v, tr ) h (1)2 (( v, tr ) ! . Ignoring higher order terms in ε (which are irrelevant for Tikhonov-Fenichelreduction) and renaming variables, one obtains the same system by settingthe second part of (cid:18) ˙ x ˙ x (cid:19) = (cid:18) P P (cid:19) · x + ε h (1)1 h (1)2 ! + · · · equal to zero, solving for x , substituting into the first part and passingto slow time. This is another proof of the fact that singular perturbationreduction and QSS reduction agree when the critical manifold is a coordinatesubspace. (The first proof was given in [12], Proposition 5.)Finally, with a view on chemical reaction networks, we address conservationlaws. Proposition 2.
Let the smooth real-valued function ψ be defined on some opensubset of U which has nonempty intersection with Φ( W ) , and assume that ψ isa first integral of system (1) for every ε . Then e ψ := ψ ◦ Φ is constant or a firstintegral of the parameterized reduced system (9) .Proof. By [16], Proposition 8 the restriction of ψ to the critical manifold Z isalso a first integral of the reduced system (4), thus Dψ ( x ) Q ( x ) h (1) ( x ) = 0 forall x ∈ Z . Therefore D e ψ ( v ) R ( v ) h (1) (Φ( v )) = Dψ (Φ( v )) D Φ( v ) R ( v ) h (1) (Φ( v ))= Dψ (Φ( v )) Q (Φ( v )) h (1) (Φ( v ))= 0 , which is the characterizing property for first integrals of system (9). The following small examples have the primary function to illustrate the argu-ments and reduction procedures from the previous subsection.8. We consider a (hypothetical) slow-fast system, with fast reaction X + X ⇋ X and slow reaction X + X ⇋ X , with associated differential equation (according to the procedure from Sub-section 3.1 below)˙ x = − k x x + k − x − εk x x + εk − x ˙ x = − k x x + k − x + 2 εk x x − εk − x ˙ x = k x x − k − x − εk x x + εk − x . The critical manifold Z is determined by Kx x = x , with K = k /k − , andwe have P = (1 , , − tr , µ ( x ) = ( − k x x + k − x ). A parameterization of Z is given by Φ : R ≥ → R , (cid:18) v v (cid:19) v v Kv v , hence to determine R ( v ) via Theorem 1(c) we have to solve R ( v ) · Kv Kv − = (cid:18) (cid:19) . By straightforward calculations one obtains R ( v ) = 11 + K ( v + v ) (cid:18) Kv − Kv − Kv Kv (cid:19) . With h (1) (Φ( v )) = ( − k Kv v + k − v ) · − the reduced system becomes (cid:18) v ′ v ′ (cid:19) = k Kv v − k − v K ( v + v ) (cid:18) − − Kv Kv (cid:19) . This is a case where the critical manifold is the graph of the rational function( x , x ) Kx x ; thus Corollary 1 would also be applicable. Moreover weget A (Φ( v )) = Dµ (Φ( v )) P (Φ( v )) = − ( k ( v + v ) + k − ) < R ≥ , hence linear stability of the critical manifold follows by Remark 3(b).9. As a non-hypothetical variant we discuss the system with the same fast re-action as in part 1, but with slow reaction X ⇋ X + X , and associated differential equation˙ x = − k x x + k − x + εk x ˙ x = − k x x + k − x ˙ x = k x x − k − x − εk x after discarding the equation for x . This is Michaelis-Menten with slowdegradation of complex to enzyme and product. Here R ( v ) is the same as inthe previous example, and h (1) (Φ( v )) = k Kv v · − . The reduced system becomes (cid:18) v ′ v ′ (cid:19) = k Kv v K ( v + v ) (cid:18) Kv − (1 + Kv ) (cid:19) , and we note a further built-in reduction: The differential equation for thereaction network admits the first integral ψ = x + x from stoichiometry,hence by Proposition 2 the reduced equation inherits the first integral e ψ = v + Kv v . Thus one ends up with a one dimensional reduced equation, asit should be.3. For contrast, consider the hypothetical slow-fast system with fast reaction2 X + 2 X ⇋ X and the same slow reaction as in part 1. The differential equation nowbecomes˙ x = − k x x + 2 k − x − εk x x + εk − x ˙ x = − k x x + 2 k − x + 2 εk x x − εk − x ˙ x = 3 k x x − k − x − εk x x + εk − x , and the critical manifold Z is given by Kx x = x , with K = k /k − , and P = (2 , , − tr . A parameterization of Z is given byΦ : (cid:18) v v (cid:19) v v Kv v . It is obvious that L = (cid:18) − (cid:19)
10s of rank two and satisfies L · P = 0. This yields, by Theorem 1(d), L · D Φ( v ) = L · v
00 3 v Kv v Kv v = (cid:18) v − v v + 4 Kv v Kv v (cid:19) . Finally the reduced system is given by v ′ = 13 v v (4 K ( v + v ) + 9 v v ) (cid:18) Kv v v − v − Kv v v (cid:19) · L · h (1) = 13 v v (4 K ( v + v ) + 9 v v ) · (cid:18) − K v v k + 12 Kv v k − − Kv v v k + 15 v k − K v v k + 27 Kv v k − Kv v k − − Kv v v k − v v k − (cid:19) . One can further reduce the dimension to one by utilizing the first integral ψ = 4 x + 5 x + 6 x from stoichiometry.In this example we could have chosen a different parameterizationΨ : (cid:18) v v (cid:19) v v Kv / v / , for Z , which directly represents Z as the graph of a function, but it may bemore convenient to work with a reduced system that has rational right handside.4. Finally we sketch an example that is motivated by mechanics, to illustratethat the range of applications does not only include reaction networks. (SeeArnold and Anosov [1] for background, and also [24].) Specifically we lookat a pair of coupled nonlinear oscillators˙ y = y + · · · ˙ y = − y + · · · ˙ y = ωy + · · · ˙ y = − ωy + · · · with irrational ω >
0, thus we are in a non-resonant scenario. Computing anormal form up to degree three and reduction by invariants y + y , y + y yields a two-dimensional system, which generically allows to decide aboutstability. But here we look at a degenerate case, with reduced equation(11) ˙ x = x ( ax + bx )˙ x = cx ( ax + bx )and parameters a < b > c <
0. The choice of signs ensures thatthe line of stationary points given by µ := ax + bx = 0 lies in the posi-tive quadrant (which is positively invariant), and also that solutions on theinvariant lines x = 0 resp. x = 0 converge to 0 as t → ∞ .11e consider (11) as fast part h (0) of a singularly perturbed system, thus wehave P ( x ) = (cid:18) x cx (cid:19) and choose Φ( v ) = (cid:18) bv − av (cid:19) , v > . A straightforward calculation shows that A (Φ( v )) = ab (1 − c ) v < v >
0, so the line µ = 0 is attracting in the positive quadrant. Moreover onemay choose L ( x ) = ( cx , − x ) and thus obtains the reduced equation(12) v ′ = − ab (1 − c ) · (cid:0) ac b (cid:1) h (1) (Φ( v ))for arbitrary small perturbation h (1) . The choice h (1) ( x ) = (cid:18) x − x (cid:19) , h (1) (Φ( v )) = (cid:18) b v − a v (cid:19) is compatible with the mechanical context under consideration here, andyields a positive stationary point for (12) as well as (11). For the originalsystem this yields the existence of an invariant torus. While the range of applications of Theorem 1 is not restricted to chemical reac-tion networks, it is natural to discuss these in greater detail: The considerationof slow and fast reactions leads to critical manifolds that consist of stationarypoints of a subnetwork, and for some relevant and familiar classes of networksexplicit parameterizations of the variety of stationary points exist. We first re-call some general facts about reaction networks and then discuss two specialclasses. The results will be illustrated by examples.
We briefly recall here the mathematical description of reaction networks accord-ing to Feinberg [5], Horn and Jackson [14] (see also the recent monograph [6] byFeinberg), and then outline the general setup for networks with fast and slowreactions, as already suggested by some illustrative examples in the previoussection.A reaction network on a set of species { X , . . . , X n } is a digraph whose nodesare finite linear combinations of species with nonnegative integer coefficients;each edge is called a reaction. Thus a node is of the form y = P ni =1 a i X i and is identified with the vector y = ( a , . . . , a n ) ∈ R n . We let x i denote theconcentration of X i and x = ( x , . . . , x n ). For each reaction (denoted by y → y ′ )we assume given a rate function w y → y ′ ( x ) ∈ R ≥ for x ∈ R n ≥ . This leads toa system of differential equations describing the evolution of the concentrationsin time:(13) ˙ x = X reactions y → y ′ w y → y ′ ( x )( y ′ − y ) , x ∈ R n ≥ . y ′ − y ∈ R n encodes the net production of each species by the occur-rence of the reaction y → y ′ . The vector subspace spanned by all the y ′ − y iscalled the stoichiometric subspace of the reaction network.It is convenient to write the system in matrix-vector form by introducingthe matrix N whose columns are the vectors y ′ − y (after fixing an order of theset of reactions). Then, with w ( x ) denoting the vector of rate functions in thesame order, (13) can be rewritten as(14) ˙ x = N w ( x ) , x ∈ R n ≥ . A frequent choice of rate function is the one from mass action kinetics , with w y → y ′ ( x ) = k y → y ′ n Y i =1 x y i i , with k y → y ′ > = 1.For the following we recall some definitions: Definition 1.
Let x, y ∈ R n and M ∈ R n × m , with columns M , . . . , M m ∈ R n .(a) For x ∈ R n> we define x y := n Y i =1 x y i i , noting that the definition may be extended to all x ∈ R n when all y i arenonnegative integers.(b) For x ∈ R n> we define x M := x M ... x M m ∈ R m , noting that the definition may be extended to all x ∈ R n when all entriesof M are nonnegative integers.(c) The Hadamard product x ◦ y is defined as the componentwise multiplicationof the two vectors x, y , i.e. x ◦ y = x · y ... x n · y n . In view of the last definition we can rewrite the reaction network for massaction kinetics in the form(15) ˙ x = N · w ( x ) = N · ( K ◦ x Y ) , K ∈ R m> ( m is the number of reactions) is a vector containing thereaction rate constants and Y ∈ R n × m is the matrix whose columns are thereactant vectors of each reaction, called the kinetic order matrix.It follows directly from (14) that any vector in the left-kernel of N defines alinear first integral, regardless of the form of w ( x ). These linear first integralsare commonly referred to as conservation laws, and their common level setsare called stoichiometric compatibility classes. If each connected component ofthe reaction network has exactly one terminal strongly connected component,then all linear first integrals of (14) arise in this way; see Feinberg and Horn[4]. Finally we recall the notion of deficiency of the reaction network, whichis defined as the number of nodes minus the rank of N minus the number ofconnected components; see e.g. Horn [15] or Feinberg [3].We turn now to a scenario with prescribed slow and fast reactions; see alsothe discussions in Heinrich and Schauer [13], Lee and Othmer [17]. The subdi-graph induced by the fast reactions is itself a reaction network with the sameset of species, which we call the fast subnetwork . We stipulate that even if somespecies are not part of any fast reaction, we still consider them as part of thefast subnetwork. We have a corresponding stoichiometric matrix N f and ratevector w f ( x ), such that, in the notation of Section 2,(16) h (0) ( x ) = N f w f ( x ) = N f · ( K f ◦ x Y f ) . Analogously, we have(17) h (1) ( x ) = N s w s ( x ) = N s · ( K s ◦ x Y s )for the slow subsystem. Keeping the notation from Section 2, we let Z be thezero set of h (0) (possibly restricted to a neighborhood e U of some point), andlet r denote the rank of Dh (0) ( x ), x ∈ Z . Then clearly rank N f ≥ r , but theinequality may be strict; see Heinrich and Schauer [13] and also Section 3 of[11]. In the present paper we will, however, restrict attention to the case whenequality holds: Blanket hypothesis.
We impose on system (16) the conditions of Proposi-tion 1(a) and the additional condition that rank N f = rank Dh (0) ( x ) = r , x ∈ Z .Due to nonnegativity of concentrations, the points of Z will be in R n ≥ . Insome instances we will require that the neighborhood e U in Proposition 1 is evena subset of R n> , and likewise we will occasionally require that the domain W of the parameterization Φ is a subset of R s> . By our assumption the zero set Z of h (0) has dimension s = n − r . Assume now that there exists a smoothparameterization Φ : W → Z with rank D Φ( v ) = s on W . 14 roposition 3. Let system (16) be given, with a parameterization Φ of thecritical manifold as in (7) , and assume the blanket hypothesis holds on Φ( W ) .Let L f ∈ R s × n be a matrix whose rows form a basis of the left-kernel of N f .Then the matrix R ( v ) in Theorem 1 is given as R ( v ) = (cid:0) L f D Φ( v ) (cid:1) − L f , and the reduced system is v ′ = (cid:0) L f D Φ( v ) (cid:1) − L f h (1) (Φ( v )) , v ∈ W. Furthermore, the column space of the matrix P ( x ) in any decomposition of h (0) ( x ) as in Proposition 1(c) equals the column space of N f .Proof. This follows from Theorem 1(d), since L f Dh (0) ( x ) = L f N f Dw f ( x ) = 0on Z .Note that by Remark 2, if rank N = rank N f , then the reduced system isgiven by v ′ = 0 . To find a function Φ( v ) which yields a parameterization of positive steady statesof the fast subnetwork, several strategies can be employed. We review here thetwo most common approaches. Throughout we use the blanket hypothesis,denoting the rank of the stoichiometric matrix N f by r , the number of species ofthe full network by n , and let s = n − r . For simplicity, we consider mass actionkinetics, although several results hold for more general classes of rate functions. As was shown in [7], non-interacting sets of species may be utilized to findrational parameterizations of the steady states, given certain conditions. Thusconsider a subset of species Y = { X i , . . . , X i r } , with the following assumptions.(i) For every fast reaction y → y ′ , both the sum of the coefficients of thespecies in Y in y and the sum of the coefficients in y ′ is at most one.This means that no pair of species in Y appear together at one side of areaction, and further no species appears with coefficient greater than 1.(ii) The rank of the submatrix of N f given by the rows i , . . . , i r is equal to r .(iii) Consider the network induced by the fast subnetwork by setting all speciesnot in Y to zero. For each species X i j in Y , there is a directed path from X i j to 0 in this induced network.In the nomenclature of [7], assumption (i) means that Y is non-interacting, (ii)means that no conservation law has support in Y , and (iii) means that there15xists a spanning tree rooted at ∗ in the appropriate digraph (see [7], Section 8for details).Let X ℓ , . . . , X ℓ s be the species not in Y . If Y satisfies (i), (ii) and (iii), thenthe components i , . . . , i r of h (0) ( x ) form a linear system in x i , . . . , x i r that hasa unique solution in terms of x ℓ , . . . , x ℓ s . Furthermore, the solution is a rationalfunction in x ℓ , . . . , x ℓ s and in the reaction rate constants k y → y ′ >
0, with allcoefficients positive [7]. The solution can be found using graphical procedures,but in practice, solving the system of linear equations is the preferred approach(see [7] for more on this).By this procedure one obtains a parameterization of the zero set Z of h (0) in s variables v i = x ℓ i , i = 1 , . . . , s . Further, clearly rank D Φ( v ) = s . In [8] someconditions are stated which guarantee that the assumptions in Proposition 1(a)are satisfied. We next consider another common scenario occurring, for instance, for so-calledcomplex balanced steady states (see Feinberg [5], Horn and Jackson [14]) andnetworks with toric steady states (see P´erez Mill´an et al. [20], M¨uller et al.[18]). In this scenario the zero set Z of h (0) in R n> agrees with the solution setof a collection of binomial equations(18) a ℓ ( k ) x u ℓ − b ℓ ( k ) x c ℓ = 0 , x ∈ R n> , ℓ = 1 , . . . , q, where u ℓ , c ℓ ∈ R n and a ℓ ( k ) , b ℓ ( k ) are polynomials in the parameters of therate functions that only attain positive values for valid k . Here, for the sakeof simplicity, we restrict attention to the case when all x i >
0. Under theseassumptions, the solution set Z to (18) equals the solution set of(19) x u ℓ − c ℓ = b ℓ ( k ) a ℓ ( k ) , x ∈ R n> , ℓ = 1 , . . . , q. The solution set of (19), if non-empty, admits a monomial parameterization ofthe following form. Let x ∗ be any fixed solution of (19) and M ∈ R n × q thematrix whose columns are u ℓ − c ℓ . If x is a solution to (19) then x M = ( x ∗ ) M . It is a classical result (see for example Lemma 3.7 in M¨uller et al. [18]) that thesolution set to this equation, and hence Z , can be parameterized in the form(20) Φ( v ) = x ∗ ◦ v B = ( x ∗ i v b i ) i =1 ,...,n , v ∈ R d> . Here d = dim ker M tr , and b , . . . , b n are the columns of a matrix B ∈ R d × n with row span equal to ker M tr and ker B tr = { } (thus d = rank B ). With aneasy computation one verifies the well-known identity(21) D Φ( v ) = diag( x ∗ ◦ v B ) B tr diag(1 /v ) , α , diag( α ) denotes the diagonal matrix with the entries of α in the diagonal, and 1 /v is defined component-wise. By this identity, therank of D Φ( v ) equals d , the rank of B , and therefore we are in the setting ofSubsection 2.2 provided that d = s . With this in place, by Proposition 3 thematrix R ( v ) in Theorem 1 becomes(22) R ( v ) = (cid:0) L f diag( x ∗ ◦ v B ) B tr diag(1 /v ) (cid:1) − L f = diag( v ) (cid:0) L f diag( x ∗ ◦ v B ) B tr (cid:1) − L f . We turn now to the special case of complex-balanced steady states for massaction kinetics. These are steady states such that for each fixed node y of thefast subnetwork, it holds that(23) X reaction y ′ → y w y ′ → y ( x )( y − y ′ ) = X reaction y → y ′ w y → y ′ ( x )( y ′ − y ) . As shown in Feinberg [5] and Horn [15], a necessary condition for complexbalanced steady states to exist is that each connected component of the fastsubnetwork is strongly connected; this is what is known as a weakly reversiblereaction network. In this case, if the parameters k y → y ′ satisfy certain algebraicconditions, then there are positive complex balanced steady states and anypositive steady state is complex balanced. Furthermore, if the deficiency of thefast subnetwork is zero, then any positive steady state is complex balanced,independent of the values of the parameters.The set Z of positive complex balanced steady states agrees with the solutionset of a collection of binomial equations of the form K ij x y i − K ji x y j = 0 , x ∈ R n> , for every pair of nodes y i , y j in the same connected component of the reactionnetwork, and with K ij and K ji positive polynomials in the parameters k y → y ′ for the reactions in the same connected component [2]. This implies that thecolumn span of the matrix M above agrees with the column span of N f , andtherefore ker M tr has rank s . As a suitable matrix B one can choose L f as inProposition 3. We thus obtain an explicit expression for the reduced system. Proposition 4.
Assume that the fast subsystem (16) has a positive complexbalanced steady state x ∗ . Then:(a) With the notation of Theorem 1 and Proposition 3, there is a parameter-ization (20) of the critical manifold with B = L f , parameter space R s> and R ( v ) = diag( v ) (cid:0) L f diag( x ∗ ◦ v L f ) L tr f (cid:1) − L f . (b) The reduced system is given by: (24) v ′ = R ( v ) · h (1) (Φ( v )) = R ( v ) · N s · (cid:0) K s ◦ ( x ∗ ) Y s ◦ v L f · Y s (cid:1) = diag( v ) (cid:0) L f diag( x ∗ ◦ v L f ) L tr f (cid:1) − L f · N s · (cid:0) K s ◦ ( x ∗ ) Y s ◦ v L f · Y s (cid:1) . roof. Part (a) is clear, while part (b) follows immediately with part (a) and h (1) (Φ( v )) = h (1) (cid:0) x ∗ ◦ v L f (cid:1) = N s · (cid:16) K s ◦ (cid:0) x ∗ ◦ v L f (cid:1) Y s (cid:17) = N s · (cid:0) K s ◦ ( x ∗ ) Y s ◦ ( v L f ) Y s (cid:1) = N s · (cid:0) K s ◦ ( x ∗ ) Y s ◦ v L f · Y s (cid:1) . In the discussion so far we restricted attention to computing a reduced system onthe critical manifold and did not address the question whether the attractivitycondition from Proposition 1(a) is satisfied. Of course, Remarks 1 and 3 areavailable, but due to our consideration of slow and fast reaction networks onealso may resort to known properties of certain classes of reaction networks.There are general attractivity results available for complex balanced mechanisms as introduced by Horn [15], which are of primary interest in our setting, since weonly require positivity of reaction constants in our considerations. By Horn [15],Theorem 4A a mechanism is complex balanced for all choices of reaction rateconstants if and only if it is weakly reversible and has deficiency zero. For thesesystems Feinberg [5] proved in Remark C.2 that every steady state is linearlyattractive (within its stoichiometric compatibility class). We therefore have:
Proposition 5.
If the fast subsystem (16) is weakly reversible and of deficiencyzero, then all non-zero eigenvalues of the Jacobian have negative real part, henceall the conditions of Proposition 1(a) hold.
We consider the following reaction network with mass action ki-netics, where the numbers k y → y ′ are written as labels of the reactions: X k −− ⇀↽ −− k X X + X k −− ⇀↽ −− k X k −− ⇀↽ −− k X + X (25) X k −−→ X X + X k −− ⇀↽ −− k X . This network is a two-component system where X , X are the unphospho-rylated and phosphorylated forms of the histidine kinase and X , X are theunphosphorylated and phosphorylated forms of the response regulator. Fur-ther we have a dead-end complex between the unphosphorylated forms of bothproteins.We now look at the slow-fast scenario where the fast reactions are those withlabels k , . . . , k , such that the fast subnetwork is X k −− ⇀↽ −− k X X + X k −− ⇀↽ −− k X k −− ⇀↽ −− k X + X X , (26) 18nd the slow reactions have labels k , k , k . With this choice of slow-fast reac-tions, we have N f = − − − − − −
10 0 1 − − , N s = − − − − and Y s = , K s = k k k . The fast reaction network in (26) has 6 nodes and three connected components,and the rank of N f is r = 3 (hence also s = 3). Therefore the deficiency is zeroand any positive steady state, that is, any element of Z , is complex balancedand a solution to a set of binomial equations. Under mass action, the steadystates of this fast subnetwork are the solutions to − k x x + k x − k x + k x = 0 , − k x x + k x = 0 , (27) k x x + k x x − k x − k x = 0 , and we easily verify that x ∗ = (cid:0) , k k , k k k k k k , , k k , (cid:1) tr is a positive steady state of the fast subnetwork. With Proposition 4(b) thereduced system can be computed. We choose(28) L f = and obtain the following parameterization of Z :Φ( v ) = x ∗ ◦ v v v L f = x ∗ ◦ v v v v v v v = v k k v k k k k k k v v k k v v v . v ′ = 1 ξ · k k k ( k k k + k k k )( − k k k k v v + k k k k v ) k ( k + k )( − k k k k v v + k k k k v ) ξk k k · ( k k k k v v − k k k k v ) , where ξ is given by ξ = k ( k k k + k k k ) k v + k ( k k k + k k k ) k v + k ( k + k )( k k k + k k k ) . In addition, we conclude by Proposition 5 that all non-zero eigenvalues of Dh (0) have negative real part on Z .Observe that the parameterization Φ( v ) is not unique. For example, choosinganother starting steady state x ∗ = (cid:16) , k k , , k k k k k k , k k k k , (cid:17) tr , we obtain the parameterizationΦ( v ) = (cid:16) v , k v k , v , k k k v k k k , k k v v k k , v (cid:17) tr , and the reduced system v ′ = − k k q ( v ) ( k k k + k k k )( k v v − k v ) v ′ = − k k q ( v ) k k ( k + k )( k v v − k v ) v ′ = k v v − k v , with q ( v ) = k k ( k k k + k k k ) v + k k k k ( k + k ) v + k ( k + k )( k k k + k k k ) . The same parameterization is obtained by eliminating x , x , x after realizingthat the set of species { X , X , X } satisfies (i), (ii) and (iii) in Subsection 3.2.1. Example 2.
If we remove the reactions with label k , k from the network(25), then the stoichiometric matrices of both the fast subnetwork and the fullnetwork have rank 3. Hence, if the reduction with a parameterized criticalmanifold is possible, the reduced system is v ′ = 0. Example 3.
We analyse the reaction network X + X k −− ⇀↽ −− k X k −−→ X k −− ⇀↽ −− k X + X k −−→ X + X (29) X k −−→ X X k −−→ X . X thekinase catalysing the phosphorylation of a substrate S with two phosphoryla-tion sites. Then X , X , X correspond to the phosphoforms with no, one,two phosphate groups respectively, and X and X are intermediate enzyme-substrate forms. Dephosphorylation proceeds without a phosphatase. We letthe fast system be all reactions involved in the conversion X ↔ X , namelythose with label k , . . . , k . Hence the reactions with label k , k are slow. Withthis choice, we have N f = − − − − − − − −
10 0 0 0 0 0 , L f = − , and h (0) ( x ) = (cid:0) − k x x − k x x + k x + k x , − k x x + k x + k x ,k x x − k x − k x , k x x + k x − k x , − k x x + k x − k x , (cid:1) tr h (1) ( x ) = (cid:0) , , , , − k x x + k x , k x x − k x (cid:1) tr . The fast network has deficiency 1, since the rank of N f is 3, and the networkhas 6 nodes and two connected components. Thus, the steady states are notcomplex balanced for all k . Instead, we observe that the set { X , X , X } satisfies (i)-(iii) in Subsection 3.2.1. Indeed, (i) and (ii) are easy to check. For(iii), the induced network obtained after setting the species not in this set tozero is 0 / / X o o (cid:15) (cid:15) X O O / / X o o and clearly there is directed path to 0 from every species. This implies that x , x , x can be solved from the system h (0) ( x ) , , = 0 to obtain the followingparameterization (where v = x , v = x , v = x )Φ : ( v , v , v ) v v k k + k v v k k k k ( k + k ) ( k v + k ) v v k k k ( k + k ) v v v , R ( v )and the reduced system using Proposition 3. We have D Φ( v ) = k v k + k k v k + k k k (2 k v + k ) v k k ( k + k ) ( k v + k ) k k v k k ( k + k ) k k v k ( k + k ) k k v k ( k + k )
00 0 1 , and using R ( v , v , v ) = (cid:0) L f D Φ( v ) (cid:1) − L f the reduced system is given by v ′ v ′ v ′ = ξ ξ k ( k k v + k k + k k ) v − (2 k k k v v + k k k v + k k k v + k k k + k k k ) ξ k ( k + k ) , where ξ = k k k v v − ( k + k ) k k v ξ = k k k v v + ( k + k ) (cid:16) k k k k v + 2 k k k k v v + k k ( k k + k k + k k ) v + k ( k + k ) k v + k ( k + k ) k (cid:17) . We next verify that the eigenvalue condition in Proposition 1(a) is satisfied.To this end we check the eigenvalues of the matrix A ( x ) = Dµ ( x ) · P ( x ) with P = − −
10 0 11 − , µ ( x ) = − k x x − k x x + k x + k x − k x x + k x + k x k x x + k x − k x . Using the Routh-Hurwitz conditions (see Gantmacher [10], Ch. V, §
6) for itscharacteristic polynomial χ A ( λ ) = λ + σ λ + σ λ + σ of degree 3, we obtain: σ = k x + k x + k x + k x + k + k + k + k σ = k k x ( x + x + x ) + ( k k + k k + k k + k k + k k ) x + k ( k + k + k ) x + k ( k + k + k ) x + ( k + k )( k + k ) + k k σ = k k k x ( x + x + x ) + k ( k k + k k + k k ) x + k k ( k + k ) x + k ( k + k )( k x + k ) . In order to verify that all eigenvalues have negative real part we use the Hurwitzconditions for polynomials of degree three: σ > , σ > , σ · σ − σ > . x with all coefficients positive, and hence are positive whenevaluated at positive values of k and x . Appendix
For convenient reference, we record a basically known result.
Lemma 1.
Let < s < n and A ∈ R n × s , B ∈ R s × n such that rank AB = s and ( AB ) = AB.
Then BA = I s .Proof. Recall that R n is the direct sum of the eigenspaces of AB for eigenvalues1 and 0. Let z , . . . , z s be a basis of the former. Then Bz , . . . , Bz s are linearlyindependent due to X λ i Bz i = 0 ⇒ X λ i ABz i = X λ i z i , hence they form a basis of R s . Finally, ABz i = z i implies ( BA ) Bz i = Bz i for1 ≤ i ≤ s . Acknowledgements.
The work of NK and SW has been supported by thebilateral project ANR-17-CE40-0036 and DFG-391322026 SYMBIONT. EF hasbeen supported by the Independent Research Fund of Denmark. SW thanks theDepartment of Mathematical Sciences of the University of Copenhagen for itshospitality during a research visit when essential parts of the present manuscriptwere devised. Likewise, EF thanks the hospitality of RWTH Aachen where thisproject was initiated.
References [1] V.I. Arnold, D.V. Anosov (eds):
Dynamical Systems I.
Springer, Berlin(1988).[2] G. Craciun, A. Dickenstein, A. Shiu, B. Sturmfels:
Toric dynamical sys-tems.
J. Symbolic Computation (11), 1551–1565 (2009).[3] M. Feinberg: Complex balancing in general kinetic systems . Arch. Ration.Mech. Anal. , 187–194 (1972).[4] M. Feinberg, F.J.M. Horn: Chemical mechanism structure and the coin-cidence of the stoichiometric and kinetic subspaces . Arch. Rational. Mech.Anal. :1, 83–97 (1977).[5] M. Feinberg: The existence and uniqueness of steady states for a classof chemical reaction networks . Arch. Ration. Mech. Anal. , 311–370(1995). 236] M. Feinberg:
Foundations of Chemical Reaction Network Theory.
Springer,Cham (2019).[7] E. Feliu, C. Wiuf:
Variable elimination in chemical reaction networks withmass action kinetics.
SIAM J. Appl. Math. , 959 –981 (2012).[8] E. Feliu, C. Lax, S. Walcher, C. Wiuf: Non-interacting species and reduc-tion of reaction networks.
In preparation.[9] N. Fenichel:
Geometric singular perturbation theory for ordinary differen-tial equations . J. Differential Equations (1), 53–98 (1979).[10] F.R. Gantmacher: Applications of the theory of matrices.
Dover, Mineola(2005).[11] A. Goeke, S. Walcher:
A constructive approach to quasi-steady state reduc-tion.
J. Math. Chem. , 2596 - 2626 (2014).[12] A. Goeke, S. Walcher, E. Zerz: Classical quasi-steady state reduction – Amathematical characterization . Physica D : 11-26 (2017).[13] R. Heinrich, M. Schauer:
Quasi-steady-state approximation in the math-ematical modeling of biochemical networks.
Math. Biosci. , 155–170(1983).[14] F. Horn, R, Jackson: General mass action kinetics.
Arch. Ration. Mech.Anal. , 81–116 (1972).[15] F. Horn: Necessary and sufficient conditions for complex balancing inchemical kinetics.
Arch. Ration. Mech. Anal. , 172–186 (1972/73).[16] C. Lax, S. Walcher: Singular perturbations and scaling.
Discrete Contin.Dyn. Syst. Ser. B (to appear). Arxiv 1807.03107 (2018).[17] C.H. Lee, H.G. Othmer:
A multi-time-scale analysis of chemical reactionnetworks: I. Deterministic systems . J. Math. Biol. , 387–450 (2009)[18] S. M¨uller, E. Feliu, G. Regensburger, C. Conradi, A. Shiu, A. Dicken-stein: Sign Conditions for Injectivity of Generalized Polynomial Maps withApplications to Chemical Reaction Networks and Real Algebraic Geometry.
Found. Comput. Math. (1), 69–97 (2016).[19] L. Noethen, S. Walcher: Tikhonov’s theorem and quasi-steady state . Dis-crete Contin. Dyn. Syst. Ser. B (3), 945–961 (2011).[20] M. Perez Millan, A. Dickenstein, A. Shiu, C. Conradi: Chemical ReactionSystems with Toric Steady States.
Bull. Math. Biology (5), 1027–1065(2012).[21] M. Stiefenhofer: Quasi-steady-state approximation for chemical reactionnetworks . J. Math. Biol. , 593–609 (1998).2422] A.N. Tikhonov: Systems of differential equations containing a small param-eter multiplying the derivative (in Russian). Math. Sb. , 575–586 (1952).[23] F. Verhulst: Methods and Applications of Singular Perturbations. BoundaryLayers and Multiple Timescale Dynamics . Springer, New York (2005).[24] S. Walcher:
On differential equations in normal form . Math. Ann.291