Tikhonov-Fenichel reduction for parameterized critical manifolds with applications to chemical reaction networks

We derive a reduction formula for singularly perturbed ordinary differential equations (in the sense of Tikhonov and Fenichel) with a known parameterization of the critical manifold. No a priori assumptions concerning separation of slow and fast variables are made, or necessary. We apply the theoretical results to chemical reaction networks with mass action kinetics admitting slow and fast reactions. For some relevant classes of such systems, there exist canonical parameterizations of the variety of stationary points; hence, the theory is applicable in a natural manner. In particular, we obtain a closed form expression for the reduced system when the fast subsystem admits complex-balanced steady states.


Introduction
A fundamental result on singular perturbation reductions, due to Tikhonov and Fenichel, allows to reduce the dimension of an ordinary differential equation with a small positive parameter in the asymptotic limit when the parameter approaches zero. This theorem has numerous applications in the sciences, in particular to chemical and biochemical reaction networks, especially quasi-steady state (QSS) for certain chemical species, and partial equilibrium approximation (PEA) for slow and fast reactions. However, the application of Tikhonov's and Fenichel's theory to reaction networks may pose some computational problems, and the purpose of the present paper is to address and resolve one of these problems.
Throughout the present work, we assume that a suitable small parameter (in a system possibly depending on several parameters) has been identified; hence, we deal with a singularly perturbed ordinary differential equation. However, we do not assume the equation to be given in separated fast and slow variables, so the usual version of the reduction theorem is not directly applicable. In applications, 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. Generally one may circumvent (and to some extent resolve) this problem by resorting to an "implicit" version of the reduction, which admits the critical submanifold of phase space as an invariant set, but this approach may also encounter computational feasibility problems.
Given this background, we derive in the present paper an explicit singular perturbation reduction that is applicable whenever a parameterization of the critical manifold is known. This reduction formalism seems particularly useful for reaction networks when the partial equilibrium approximation is applicable, since in many instances varieties of stationary points admit a canonical parameterization. We study this scenario in detail, and our results include a general closed form reduction formula, based on stoichiometry alone, for fast subsystems that are complex balanced, as well as a discussion of linear attractivity properties of slow manifolds.
The paper is organized as follows. After some preliminary work (mostly recalling notions and results from the literature) we derive in Sect. 2 a general formula for Tikhonov-Fenichel reduction when a (possibly local) parameterization of the critical manifold is given, and also consider some special cases. We illustrate the procedure by some small examples, briefly indicating that the range of applications is not restricted to chemical reaction networks. The setting of reaction networks is studied in Sect. 3. To finish the paper we discuss some special examples.

A Reduction Formula
We consider the singular perturbation reduction of ordinary differential equations, with no a priori assumption on separated (slow and fast) variables. Thus let U ⊆ R n be open, ε 0 > 0, and let h be a smooth function in some neighborhood of U × [0, ε 0 ), with values in R n . This defines a parameter-dependent system of ordinary differential equations, viz.

Journal of Nonlinear Sciencė
x = h (0) (x) + εh (1) and rewritten in slow time scale τ = εt we have a singularly perturbed system Here h (0) is called the fast part and h (1) the slow part of either system. We focus on the behavior of (1), (2) as ε → 0, and we will restrict attention to scenarios for which, modulo a coordinate transformation, the classical singular perturbation theorems of Tikhonov (1952) and Fenichel (1979) are applicable. (To be specific, we refer to the version of Tikhonov's theorem as given in Verhulst 2005, Theorem 8.1 ff.; see also Goeke and Walcher 2014, Sect. 2.1.) The focus of the present paper is on singularly perturbed systems which are not in "Tikhonov standard form" with separated slow and fast variables, but admit a transformation to such a form. Finding reductions explicitly is relevant for various applications. In Sect. 3 we will discuss in detail the application to chemical reaction networks.
The earliest result about reduction in a coordinate-free setting is due to Fenichel (see Fenichel 1979, Sect. 5, in particular Lemma 5.4), who discussed the case when the critical manifold is given explicitly as the graph of some function. Generally there exist intrinsic conditions for the existence of a coordinate transformation to Tikhonov standard form (see Fenichel 1979, Sect. 5 andNoethen andWalcher 2011, Proposition 2.1), but such a transformation frequently cannot be obtained in explicit form. A general implicit reduction procedure was developed in Goeke and Walcher (2014); Noethen and Walcher (2011). Building on this, in the present paper we derive a version of the reduced system that can be computed explicitly when a parameterization of the critical manifold is known.

Review: Tikhonov-Fenichel Reduction
We recall some basic results on coordinate-independent Tikhonov-Fenichel reduction from Goeke and Walcher (2014); Noethen and Walcher (2011) which are applicable whenever the critical manifold is given implicitly. In particular we refer to Goeke and Walcher (2014), Theorem 1 and the subsequent remarks. Proposition 1 Let system (1) be given, and denote by V(h (0) ) the zero set of h (0) . Moreover, let 0 < r < n and set s := n − r > 0.
(a) Assume that a ∈ V(h (0) ) has the following properties.
• There exists a neighborhood U of a such that rank Dh (0) • For all x ∈ Z the nonzero eigenvalues of Dh (0) (x) have real part < 0.
Then in some neighborhood of a there exists an invertible coordinate transformation from (1) to Tikhonov standard forṁ with separated slow and fast variables; moreover, the fast system satisfies a linear stability condition.
moreover rank P(a) = r, rank Dμ(a) = r and Here the entries of μ may be taken as any r entries of h (0) that are functionally independent at a. (d) The following system (in slow time) is defined on U , and admits Z as an invariant set: with The restriction of this system to Z corresponds to the reduced equation in Tikhonov's theorem.
We refer to Z as the local critical manifold (or local asymptotic slow manifold) of system (1). The decomposition (3) for smooth vector fields is a consequence of the implicit function theorem, and thus in general an explicit determination may not be possible. But the decomposition can be determined algorithmically whenever the right-hand side of (1) is polynomial or rational; see Goeke and Walcher (2014).

Remark 1 (a) Note that
has full rank on Z , Dh (0) (x) and P(x) have the same column space. (b) The eigenvalues of A(x), x ∈ Z , are the nonzero eigenvalues of Dh (0) (x) whenever the latter has rank r ; see Goeke and Walcher (2014), Remark 3.
(c) We call the projection operator of the reduction. For each x this is a linear projection of rank s = n − r which sends every element of R n to its kernel component from the kernel-image decomposition with respect to Dh (0) (x). (In the special case when the critical manifold is the graph of some function, a projection matrix was determined in Fenichel 1979, Lemma 5.4.) (d) Formally system (4) is defined whenever A(x) is invertible, and by Fenichel's results it corresponds to a reduced system as ε → 0 whenever all eigenvalues of A(x) have nonzero real part (normal hyperbolicity).
Remark 2 The reduced system may just have the form x = 0; in particular this occurs in the following scenario: h (0) always admits n − s independent first integrals near any point of Z , and locally every point of Z is uniquely determined as an intersection of Z with suitable level sets of these first integrals; see Goeke and Walcher (2014)), Sect. 2.3. Now, in the special case when h (1) admits the same first integrals, then Z as well as every intersection of Z with level sets is invariant for the reduced equation, meaning that every point of Z is invariant, thus stationary. However, the only information to be gained from x = 0 for small ε > 0 is that system (2) restricted to the invariant manifold has right-hand side of order ε 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 singular perturbation reduction, the critical manifold Z is given only implicitly via the zeros 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 is natural to search for simplified reduction procedures in special circumstances. One notable scenario appears when a parameterization for the critical manifold is explicitly known, and we will next discuss reduction in this case.

Parameterized Critical Manifolds
We keep the assumptions and notation from Proposition 1, in particular the decomposition (3), the s-dimensional local critical manifold Z (being the zero set of μ, as well as of h (0) ), and the reduced system Now assume that there is an open set W ⊆ R s and a smooth parameterization Then every solution x(t) of (6) with initial value in (W ) can be written in the form for t in some neighborhood of 0, and differentiation yields 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 (c) The matrix R(v) is uniquely determined by the conditions and therefore can be obtained from the matrix equation For every x ∈ Z let L(x) ∈ R s×n be of full rank s and such that L(x)Dh (0)

and the reduced system, in parameterized form, is given by
Proof For every v ∈ W one has h (0) ( (v)) = 0, and by differentiation Thus the image of D (v) is contained in the kernel of Dh (0) ( (v)), and these two 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, and the latter has full rank. Therefore, every column of Q( (v)) is a unique linear combination of the columns of D (v). Rewritten in matrix language, this is the assertion of part (a). Part (b) is now obvious from Eq. (8) and injectivity of 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 Eq. (5)), and the fact that D (v) is an injective linear map. The second condition follows from the fact that is a projection of rank s, by using Lemma 1 in "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)) span the image (see Goeke and Walcher 2014 for more details).
To prove (d), first recall from Remark 1 that Using now the second condition in part (c), we have which leads to the asserted expression.

Remark 3 (a) The second equation in part (c) of the theorem shows that R(v) is (according to definition) a left inverse of D (v).
This observation allows to obtain (9) directly from (8). As a consequence, one obtains the identities but we note that in general R(v) will not be the Moore-Penrose inverse of D (v). (b) The characterization in part (d) of the theorem shows that R(v) can be computed by standard linear algebra: One only needs to determine the left kernel of Dh (0) (x) to find L(x), and then compute products and inverses of certain matrices. There is no need for explicit knowledge of the projection matrix Q, or of the matrix P from the decomposition, as part (c) might suggest. However, the column space of On the other hand, knowledge of P and μ seems indispensable for the computation of A(x) = Dμ(x)P(x), and of A( (v)). Note that the eigenvalues of the latter provide direct information on the stability of the critical manifold; see Remark 1 (b).
We consider some special cases in more detail.

Corollary 1 (a) Assume that
and partition In these expressions the argument of D 1 and D 2 is v and the argument of P 1 and P 2 is (v). (b) (Fenichel 1979, Lemma 5.4:) In the special case when 1 (v) = v, thus the critical manifold is the graph of 2 , we get Proof With R i given as above one verifies R 1 D 1 + R 2 D 2 = I s and R 1 P 1 + R 2 P 2 = 0 by direct computation. Rewriting, one obtains which 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 condition on the derivative. Moreover, by local invertibility of 1 there exists a parameterization In the yet more special case that 1 (v) = v and 2 (v) = 0 the procedure yields the familiar quasi-steady-state reduction. Indeed, in this case one has μ(x) = x 2 for x = (x 1 , x 2 ) tr and x 1 ∈ R s , thus R 1 = I s , R 2 = −P 1 P −1 2 and the reduced equation is .
Ignoring higher order terms in ε (which are irrelevant for Tikhonov-Fenichel reduction) and renaming variables, one obtains the same system by setting the second part of equal to zero, solving for x 2 , substituting into the first part and passing to slow time. This is another proof of the fact that singular perturbation reduction and QSS reduction agree when the critical manifold is a coordinate subspace. (The first proof was given in Goeke et al. 2017, Proposition 5.) Finally, with a view on chemical reaction networks, we address conservation laws.

Proposition 2 Let the smooth real-valued function ψ be defined on some open subset of U which has non-empty intersection with (W )
, and assume that ψ is a first integral of system (1) for every ε. Then ψ := ψ • is constant or a first integral of the parameterized reduced system (9).
Proof By Lax and Walcher (2018), Proposition 8 the restriction of ψ to the critical manifold Z is also a first integral of the reduced system (4) which is the characterizing property for first integrals of system (9).

Illustrative Examples
The following small examples have the primary function to illustrate the arguments and reduction procedures from the previous subsection.
Example 1 We consider a (hypothetical) slow-fast system, with fast reaction and slow reaction with associated differential equation (according to the procedure from Sect. 3.1) By straightforward calculations, one obtains This is a case where the critical manifold is the graph of the rational function (x 1 , x 2 ) → K x 1 x 2 ; thus Corollary 1 would also be applicable. Moreover, we get hence linear stability of the critical manifold follows by Remark 3(b).
Example 2 As a non-hypothetical variant, we discuss the system with the same fast reaction as in Example 1, but with slow reaction and associated differential equatioṅ after discarding the equation for x 4 . This is Michaelis-Menten with slow degradation of complex to enzyme and product. Here R(v) is the same as in the previous example, and and we note a further built-in reduction: The differential equation for the reaction network admits the first integral ψ = x 1 +x 3 from stoichiometry; hence, by Proposition 2 the reduced equation inherits the first integral ψ = v 1 + K v 1 v 2 . Thus one ends up with a one-dimensional reduced equation, as it should be.
Example 3 For contrast, consider the hypothetical slow-fast system with fast reaction 2X 1 + 2X 2 3X 3 and the same slow reaction as in Example 1. The differential equation now becomeṡ and the critical manifold Z is given by K x 2 1 x 2 2 = x 3 3 , with K = k 1 /k −1 , and P = (2, 2, −3) tr . A parameterization of Z is given by It is obvious that is of rank two and satisfies L · P = 0. This yields, by Theorem 1(d), Finally the reduced system is given by One can further reduce the dimension to one by utilizing the first integral ψ = 4x 1 + 5x 2 + 6x 3 from stoichiometry.
In this example we could have used Remark 4(a) and chosen a different parameterization which directly represents Z as the graph of a function. But it seems more natural and convenient to work with a reduced system that has rational right-hand side. Moreover, the second parameterization may obscure the fact that Z is an algebraic variety.
Example 4 Finally we sketch an example that is motivated by mechanics, to indicate that the range of applications does not only include reaction networks. (See Arnold and Anosov 1988 for background, and also Walcher 1991.) Specifically we look at a pair of coupled nonlinear oscillators ⎛ with irrational ω > 0; thus, we are in a non-resonant scenario, f 3 a homogeneous polynomial vector field of degree three, and "t.h.o." denoting terms of higher order. We assume that the system is in normal form up to degree three, with the special form where a < 0, b > 0, c < 0 are real constants and p 1 , p 2 are linear combinations of y 2 1 + y 2 2 and y 2 3 + y 2 4 . Reduction of the degree three Taylor polynomial by the invariants y 2 1 + y 2 2 , y 2 3 + y 2 4 yields a two-dimensional system, given bẏ The reduced system is degenerate in our case, admitting a line of stationary points given by μ := ax 1 + bx 2 = 0. The choice of signs ensures that this line lies in the positive quadrant (which is positively invariant), and also that solutions on the invariant lines x 1 = 0 resp. x 2 = 0 converge to 0 as t → ∞.
We consider (11) as fast part h (0) of a singularly perturbed system, thus we have for arbitrary small perturbation h (1) . The choice is compatible with the mechanical context under consideration here and yields a positive stationary point for (12) as well as (11). For the original system, this yields the existence of an invariant torus.

Applications to Reaction Networks
While the range of applications of Theorem 1 is not restricted to chemical reaction networks, it is natural to discuss reduction of these reaction networks, for the following reason: The consideration of slow and fast reactions leads to critical manifolds that consist of stationary points of a subnetwork. Extensive previous work on parametrizations of the positive part of the variety of stationary points in chemical reaction networks (see, e.g., Horn 1972;Craciun et al. 2009;Müller et al. 2016) shows that parametrizations exist and can be explicitly found for several relevant and familiar classes of networks. Thus the results of the previous section can often be applied and reduced systems can be explicitly obtained. In particular, in the case where the fast subnetwork admits a complex-balanced steady state (see below), a closed formula for the reduced system can be readily expressed in terms of the stoichiometry of the network alone (see Proposition 4). We first recall some general facts about reaction networks and then discuss two special classes. The results will be illustrated by examples.

Reaction Networks
We briefly recall here the mathematical description of reaction networks according to Feinberg (1995), Horn and Jackson (1972) (see also the recent monograph by Feinberg (2019)), and then outline the general setup for networks with fast and slow reactions, as already suggested by some illustrative examples in the previous section.
A reaction network on a set of species {X 1 , . . . , X n } is a digraph whose nodes are finite linear combinations of species with nonnegative integer coefficients; each edge is called a reaction. Thus a node is of the form y = n i=1 a i X i and is identified with the vector y = (a 1 , . . . , a n ) ∈ R n . We let x i denote the concentration of X i and x = (x 1 , . . . , x n ). For each reaction (denoted by y → y ) we assume given a rate function w y→y (x) ∈ R ≥0 for x ∈ R n ≥0 . This leads to a system of differential equations describing the evolution of the concentrations in time: Note that y − y ∈ R n encodes the net production of each species by the occurrence of the reaction y → y . The vector subspace spanned by all the y − y is called the stoichiometric subspace of the reaction network. It is convenient to write the system in matrix-vector form by introducing the matrix N whose columns are the vectors y − y (after fixing an order of the set of reactions). Then, with w(x) denoting the vector of rate functions in the same order, (13) can be rewritten asẋ A frequent choice of rate function is the one from mass action kinetics, with with k y→y > 0 called reaction rate constants and using the convention 0 0 = 1. For the following we recall some definitions: Definition 1 Let x, y ∈ R n and M ∈ R n×m , with columns M 1 , . . . , M m ∈ R n .
(a) For x ∈ R n >0 we define noting that the definition may be extended to all x ∈ R n when all y i are nonnegative integers.
noting that the definition may be extended to all x ∈ R n when all entries of M are nonnegative integers. (c) The Hadamard product x • y is defined as the component-wise multiplication of the two vectors x, y, i.e., In view of the last definition we can rewrite the reaction network for mass action kinetics in the formẋ where K ∈ R m >0 (m is the number of reactions) is a vector containing the reaction rate constants and Y ∈ R n×m is the matrix whose columns are the reactant vectors of each reaction, called the kinetic order matrix.
It follows directly from (14) that any vector in the left kernel of N defines a linear first integral, regardless of the form of w(x). These linear first integrals are commonly referred to as conservation laws, and their common level sets are called stoichiometric compatibility classes. If each connected component of the reaction network has exactly one terminal strongly connected component, then all linear first integrals of (14) arise in this way; see Feinberg and Horn (1977). Finally we recall the notion of deficiency of the reaction network, which is defined as the number of nodes minus the rank of N minus the number of connected components; see, e.g., Horn (1972) or Feinberg (1972. We turn now to a scenario with prescribed slow and fast reactions; see also the discussions in Heinrich and Schauer (1983), Lee and Othmer (2009). The subdigraph induced by the fast reactions is itself a reaction network with the same set of species, which we call the fast subnetwork. We stipulate that even if some species are not part of any fast reaction, we still consider them as part of the fast subnetwork. We have a corresponding stoichiometric matrix N f and rate vector w f (x), such that, in the notation of Sect. 2, Analogously, we have for the slow subsystem. Keeping the notation from Sect. 2, we let Z be the zero set of h (0) (possibly restricted to a neighborhood U of some point), and let r denote the rank of Dh (0) (x), x ∈ Z . Then clearly rank N f ≥ r , but the inequality may be strict; see Heinrich and Schauer (1983) and also Sect. 3 of Goeke and Walcher (2014). In the present paper we will, however, restrict attention to the case when equality holds: Blanket hypothesis We impose on system (16) the conditions of Proposition 1(a) and the additional condition that rank N f = rank Dh (0) Due to nonnegativity of concentrations, the points of Z will be in R n ≥0 . In some instances we will require that the neighborhood U in Proposition 1 is even a subset of R n >0 , and likewise we will occasionally require that the domain W of the parameterization is a subset of R s >0 . By our assumption the zero set Z of h (0) has dimension s = n − r . Assume now that there exists a smooth parameterization

Proposition 3 Let system (16) be given, with a parameterization
of the critical 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 Note that by Remark 2, if rank N = rank N f , then the reduced system is given by v = 0.

Canonical Parameterizations for Some Classes
To find a function (v) which yields a parameterization of positive steady states of the fast subnetwork, several strategies can be employed. We review here the two most common approaches. Throughout we use the blanket hypothesis, denoting the rank of the stoichiometric matrix N f by r , the number of species of the full network by n, and let s = n − r . For simplicity, we consider mass action kinetics, although several results hold for more general classes of rate functions.

Non-interacting Sets and Rational Parameterizations
As was shown in Feliu and Wiuf (2012), non-interacting sets of species may be utilized to find rational parameterizations of the steady states, given certain conditions. Thus consider a subset of species Y = {X i 1 , . . . , X i r }, with the following assumptions.
(i) For every fast reaction y → y , both the sum of the coefficients of the species in Y in y and the sum of the coefficients in y are at most one. This means that no pair of species in Y appear together at one side of a reaction, and further no species appears with coefficient greater than 1. (ii) The rank of the submatrix of N f given by the rows i 1 , . . . , i r is equal to r . (iii) Consider the network induced by the fast subnetwork by setting all species not 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 Feliu and Wiuf (2012), assumption (i) means that Y is noninteracting, (ii) means that no conservation law has support in Y, and (iii) means that there exists a spanning tree rooted at * in the appropriate digraph (see Feliu and Wiuf 2012, Sect. 8 for details). Let X 1 , . . . , X s be the species not in Y. If Y satisfies (i), (ii) and (iii), then the components i 1 , . . . , i r of h (0) (x) form a linear system in x i 1 , . . . , x i r that has a unique solution in terms of x 1 , . . . , x s . Furthermore, the solution is a rational function in x 1 , . . . , x s and in the reaction rate constants k y→y > 0, with all coefficients positive (Feliu and Wiuf 2012). The solution can be found using graphical procedures, but in practice, solving the system of linear equations is the preferred approach (see Feliu and Wiuf 2012 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 Feliu et al. (2019) some conditions are stated which guarantee that the assumptions in Proposition 1(a) are satisfied.

Monomial Parameterizations and Deficiency Zero Networks
We next consider another common scenario occurring, for instance, for so-called complex-balanced steady states (see Feinberg 1995;Horn and Jackson 1972) and networks with toric steady states (see Perez Millan et al. 2012;Müller et al. 2016). In this scenario the zero set Z of h (0) in R n >0 agrees with the solution set of a collection of binomial equations where u , c ∈ R n and a (k), b (k) are polynomials in the parameters of the rate functions that only attain positive values for valid k. Here, for the sake of simplicity, we restrict attention to the case when all x i > 0. Under these assumptions, the solution set Z to (18) equals the solution set of The solution set of (19), if non-empty, admits a monomial parameterization of the following form. Let x * be any fixed solution of (19) and M ∈ R n×q the matrix whose columns are u − c . If x is a solution to (19), then It is a classical result (see for example Lemma 3.7 in Müller et al. 2016) that the solution set to this equation, and hence Z , can be parameterized in the form Here d = dim ker M tr , and b 1 , . . . , b n are the columns of a matrix B ∈ R d×n with row span equal to ker M tr and ker B tr = {0} (thus d = rank B). With an easy computation one verifies the well-known identity where, for a vector α, diag(α) denotes the diagonal matrix with the entries of α in the diagonal, and 1/v is defined component-wise. By this identity, the rank of D (v) equals d, the rank of B, and therefore we are in the setting of Sect. 2.2 provided that d = s. With this in place, by Proposition 3 the matrix R(v) in Theorem 1 becomes We turn now to the special case of complex-balanced steady states for mass action kinetics. These are steady states such that for each fixed node y of the fast subnetwork, it holds that reaction y →y w y →y (x)(y − y ) = reaction y→y w y→y (x)(y − y).
As shown in Feinberg (1995) and Horn (1972), a necessary condition for complexbalanced steady states to exist is that each connected component of the fast subnetwork is strongly connected; this is what is known as a weakly reversible reaction network. In this case, if the parameters k y→y satisfy certain algebraic conditions, then there are positive complex-balanced steady states and any positive steady state is complex balanced. Furthermore, if the deficiency of the fast 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 solution set of a collection of binomial equations of the form for every pair of nodes y i , y j in the same connected component of the reaction network, and with K i j and K ji positive polynomials in the parameters k y→y for the reactions in the same connected component (Craciun et al. 2009). This implies that the column span of the matrix M above agrees with the column span of N f , and therefore ker M tr has rank s. As a suitable matrix B one can choose L f as in Proposition 3. We thus obtain an explicit expression for the reduced system.

Proposition 4
Assume that the fast subsystem (16) has a positive complex-balanced steady state x * . Then: (a) With the notation of Theorem 1 and Proposition 3, there is a parameterization (20) of the critical manifold with B = L f , parameter space R s >0 and The reduced system is given by: Proof Part (a) is clear, while part (b) follows immediately with part (a) and

Attractivity of the Critical Manifold
In the discussion so far we restricted attention to computing a reduced system on the critical manifold and did not address the question whether the attractivity condition from Proposition 1(a) is satisfied. Of course, Remarks 1 and 3 are available, but due to our consideration of slow and fast reaction networks one also 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 (1972), which are of primary interest in our setting, since we only require positivity of reaction constants in our considerations. By Horn (1972), Theorem 4A a mechanism is complex balanced for all choices of reaction rate constants if and only if it is weakly reversible and has deficiency zero. For these systems Feinberg (1995) proved in Remark C.2 that every steady state is linearly attractive (within its stoichiometric compatibility class). We therefore have: Proposition 5 If the fast subsystem (16) is weakly reversible and of deficiency zero, then all nonzero eigenvalues of the Jacobian have negative real part, hence all the conditions of Proposition 1(a) hold.

Examples
Example 5 We consider the following reaction network with mass action kinetics, where the numbers k y→y are written as labels of the reactions: This network is a two-component system where X 1 , X 2 are the unphosphorylated and phosphorylated forms of the histidine kinase and X 3 , X 4 are the unphosphorylated and phosphorylated forms of the response regulator. Further we have a dead-end complex between the unphosphorylated forms of both proteins. We now look at the slow-fast scenario where the fast reactions are those with labels k 1 , . . . , k 6 , such that the fast subnetwork is and the slow reactions have labels k 7 , k 8 , k 9 . With this choice of slow-fast reactions, we have 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 zero and any positive steady state, that is, any element of Z , is complex balanced and a solution to a set of binomial equations. Under mass action, the steady states of this fast subnetwork are the solutions to and we easily verify that and obtain the following parameterization of Z : Using Eq. (24), the reduced system is found: X 2 , X 5 , X 6 correspond to the phosphoforms with no, one, two phosphate groups, respectively, and X 3 and X 4 are intermediate enzyme-substrate forms. Dephosphorylation proceeds without a phosphatase. We let the fast system be all reactions involved in the conversion X 2 ↔ X 5 , namely those with label k 1 , . . . , k 6 . Hence the reactions with label k 7 , k 8 are slow. With this choice, we have and h (0) (x) = − k 1 x 1 x 2 − k 5 x 1 x 5 + k 2 x 3 + k 4 x 4 , −k 1 x 1 x 2 + k 2 x 3 + k 6 x 5 , − k 5 x 1 x 5 + k 4 x 4 − k 6 x 5 , 0 tr h (1) (x) = 0, 0, 0, 0, −k 7 x 1 x 5 + k 8 x 6 , k 7 x 1 x 5 − k 8 x 6 tr .
The three expressions on the left side of each inequality are polynomials in the parameters and x with all coefficients positive, and hence are positive when evaluated at positive values of k and x.