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

TL;DR

This paper develops a Tikhonov-Fenichel reduction method for singularly perturbed ODEs with known critical manifolds, applicable to chemical reaction networks, especially when the fast subsystem has complex balanced steady states.

Contribution

It introduces a reduction formula that does not require separation of variables and applies to chemical networks with known stationary point parameterizations.

Findings

Derived a reduction formula for singularly perturbed ODEs.

Applied the theory to chemical reaction networks with mass action kinetics.

Obtained a closed form for the reduced system in specific cases.

Abstract

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.

Equations

$$\dot{x}=h^{(0)}(x)+\epsilon h^{(1)}(x)+{\epsilon }^2\dots ,x\in U,\epsilon \ge 0,$$

$$x^{\prime }=\frac{1}{\epsilon }{h}^{(0)}(x)+h^{(1)}(x)+\epsilon \dots ,x\in U,\epsilon \ge 0.$$

$${\mathbb{R}}^n=\mathrm{Ker}D{h}^{(0)}(x)\oplus \mathrm{Im}D{h}^{(0)}(x).$$

$$\begin{array}[]{rcl}\dot{y}_{1}&=&\varepsilon f_{1}(y_{1},y_{2})+O(\varepsilon^{2})\\ \dot{y}_{2}&=&\phantom{\varepsilon}f_{2}(y_{1},y_{2})+O(\varepsilon)\\ \end{array}$$

$$h^{(0)}(x)=P(x)\mu (x),\text{for all }x\in Z;$$

$$Z=\mathcal{V}(\mu )\cap \tilde{U}.$$

$$x^{\prime }=\left(I_n-P(x)A(x{)}^{-1}D\mu (x)\right)h^{(1)}(x),$$

$$A(x):=D\mu (x)P(x).$$

$$Q(x):=I_n-P(x)A(x{)}^{-1}D\mu (x)$$

$$x^{\prime }=Q(x)h^{(1)}(x)\text{ on }Z.$$

$$\Phi :W\to Z,\mathrm{rank}D\Phi (v)=s\text{ for all }v\in W.$$

$$x(t)=\Phi (v(t)),$$

$$D\Phi (v(t))v^{\prime }(t)=x^{\prime }(t)=Q(\Phi (v(t)))\cdot h^{(1)}(\Phi (v(t))).$$

$$Q(\Phi (v))=D\Phi (v)\cdot R(v).$$

$$v^{\prime }=R(v)\cdot h^{(1)}(\Phi (v)).$$

$$R(v)\cdot P(\Phi (v))=0\text{and }R(v)\cdot D\Phi (v)=I_s,$$

$$R(v)\cdot \left(D\Phi (v)|P(\Phi (v))\right)=\left(I_s|0\right)$$

$$R(v)=\big{(}L^{*}(v)\,D\Phi(v)\big{)}^{-1}L^{*}(v),$$

$$v^{\prime}=\big{(}L^{*}(v)\,D\Phi(v)\big{)}^{-1}L^{*}(v)\,h^{(1)}(\Phi(v)).$$

$$D{h}^{(0)}(\Phi (v))D\Phi (v)=0.$$

$$\Lambda (v)L(\Phi (v))D\Phi (v)=I_s,$$

$$\Lambda(v)=\big{(}L(\Phi(v))D\Phi(v)\big{)}^{-1},$$

$$\Phi (v)=\left(\begin{array}{c}{\Phi }_1(v)\\ {\Phi }_2(v)\end{array}\right),\text{ with }{\Phi }_1(v)\in {\mathbb{R}}^s\text{ and }D{\Phi }_1(v)\text{ invertible, for all }v\in W,$$

$$P(x)=\left(\begin{array}{c}{P}_1(x)\\ P_2(x)\end{array}\right)\text{ with }{P}_1(x)\in {\mathbb{R}}^{s\times r}.$$

$$R(v)=\left(R_1(v)|R_2(v)\right)$$

$$\begin{array}[]{rcl}R_{1}(v)&=&I_{s}-P_{1}\left(D\Phi_{2}D\Phi_{1}^{-1}P_{1}-P_{2}\right)^{-1}D\Phi_{2}D\Phi_{1}^{-1}\\ R_{2}(v)&=&P_{1}\left(D\Phi_{2}D\Phi_{1}^{-1}P_{1}-P_{2}\right)^{-1}.\end{array}$$

$$\begin{array}[]{rcl}R_{1}(v)&=&I_{s}-P_{1}\left(D\Phi_{2}P_{1}-P_{2}\right)^{-1}D\Phi_{2}\\ R_{2}(v)&=&P_{1}\left(D\Phi_{2}P_{1}-P_{2}\right)^{-1}.\end{array}$$

$$R_{1}D{\Phi }_1+R_{2}D{\Phi }_2=I_s\text{ and }{R}_1{P}_1+R_2{P}_2=0$$

$$\left(R_1|R_2\right)\left(\begin{array}{cc}D{\Phi }_1& P_1\\ D{\Phi }_2& P_2\end{array}\right)=\left(I_s|0\right),$$

$$v^{\prime }=\left(I_s|-P_1((v,0{)}^{\mathrm{tr}})P_2^{-1}((v,0{)}^{\mathrm{tr}})\right)\left(\begin{array}{c}{h}_1^{(1)}((v,0{)}^{\mathrm{tr}})\\ h_2^{(1)}((v,0{)}^{\mathrm{tr}})\end{array}\right).$$

Full text

MSC (2010): 92C45, 34E15, 80A30, 13P10

Key words: Singular perturbation, critical manifold, chemical reaction network, complex balancing

1 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.

The paper is organized as follows. After some preliminary work (mostly recalling notions and results from the literature) we derive in Section 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. However, 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. This setting is discussed in detail in Section 3, and our results include a closed form reduction formula for fast subsystems that are complex balanced, as well as a discussion of linear attractivity properties of slow manifolds. To finish the paper we discuss some examples.

2 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\subseteq\mathbb{R}^{n}$ be open, $\varepsilon_{0}>0$, and let $h$ be a smooth function in some neighborhood of $U\times[0,\varepsilon_{0})$, with values in $\mathbb{R}^{n}$. This defines a parameter-dependent system of ordinary differential equations, viz.

[TABLE]

and rewritten in slow time scale $\tau=\varepsilon t$ we have a singularly perturbed system

[TABLE]

Here $h^{(0)}$ is called the fast part and $h^{(1)}$ the slow part of either system. We focus on the behavior of the (1), (2) as $\varepsilon\to 0$, and we will restrict attention to scenarios for which, modulo a coordinate transformation, the classical singular perturbation theorems of Tikhonov [22] and Fenichel [9] are applicable. (To be specific, 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 conditions for the existence of a coordinate transformation to “Tikhonov standard form” with separated slow and fast variables, such a transformation cannot generally be obtained in explicit form. We first recall a general implicit reduction procedure developed in [11, 19], and then present, as a new result, a version of the reduced system that can be computed explicitly when a parameterization of the critical manifold is known.

2.1 Review: Tikhonov–Fenichel reduction

We recall the essential results on coordinate-independent Tikhonov-Fenichel reduction from [11, 19]; in particular we refer to [11], Theorem 1 and the subsequent remarks.

Proposition 1**.**

Let system (1) be given, and denote by ${\mathcal{V}}(h^{(0)})$ the zero set of $h^{(0)}$. Moreover let $0<r<n$ and set $s:=n-r>0$.

1. (a)

Assume that $a\in{\mathcal{V}}(h^{(0)})$ has the following properties.

• •

There exists a neighborhood $\widetilde{U}$ of $a$ such that ${\rm rank}\,Dh^{(0)}(x)=r$ for all $x\in Z:={\mathcal{V}}(h^{(0)})\cap\widetilde{U}$; in particular $Z$ is an $s$-dimensional submanifold of $\mathbb{R}^{n}$.

• •

For all $x\in Z$ there is a direct sum decomposition

[TABLE]

• •

For all $x\in 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 form

[TABLE]

with separated slow and fast variables; moreover the fast system satisfies a linear stability condition. 2. (b)

Conversely, the conditions in part (a) are necessary for the existence of a local coordinate transformation to Tikhonov standard form. 3. (c)

One may choose $\widetilde{U}$ such that there exists a product decomposition with functions $\mu(x)$ taking values in $\mathbb{R}^{r\times 1}$, $P(x)$ taking values in $\mathbb{R}^{n\times r}$, such that

[TABLE]

moreover ${\rm rank}\ P(a)=r$, ${\rm rank}\ D\mu(a)=r$ and

[TABLE]

Here the entries of $\mu$ may be taken as any $r$ entries of $h^{(0)}$ that are functionally independent at $a$. 4. (d)

The following system (in slow time) is defined on $\widetilde{U}$, and admits $Z$ as an invariant set:

[TABLE]

with

[TABLE]

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).

Remark 1**.**

1. (a)

Note that $Dh^{(0)}(x)=P(x)D\mu(x)$ on $Z$, due to $\mu(x)=0$. Since $P(x)$ has full rank on $Z$, $Dh^{(0)}(x)$ and $P(x)$ have the same column space. 2. (b)

The eigenvalues of $A(x)$, $x\in Z$, are the nonzero eigenvalues of $Dh^{(0)}(x)$ whenever the latter has rank $r$; see [11], Remark 3. 3. (c)

We call

[TABLE]

the projection operator of the reduction. For each $x$ this is a linear projection of rank $s=n-r$ which sends every element of $\mathbb{R}^{n}$ to its kernel component from the kernel-image decomposition with respect to $Dh^{(0)}(x)$. 4. (d)

Formally system (4) is defined whenever $A(x)$ is invertible, and by Fenichel’s results it corresponds to a reduced system as $\varepsilon\to 0$ whenever all eigenvalues of $A(x)$ have nonzero real part (normal hyperbolicity).

Remark 2**.**

The reduced system may just have the form $x^{\prime}=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 [11], Subsection 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^{\prime}=0$ for small $\varepsilon>0$ is that system (2) restricted to the invariant manifold has right hand side of order $\varepsilon$ or higher. (Generally the reduced system (4) in slow time represents only the $O(1)$ term in $\varepsilon$.)

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 $\mathbb{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.

2.2 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 $\mu$, as well as of $h^{(0)}$), and the reduced system

[TABLE]

Now assume that there is an open set $W\subseteq\mathbb{R}^{s}$ and a smooth parameterization

[TABLE]

Then every solution $x(t)$ of (6) with initial value in $\Phi(W)$ can be written in the form

[TABLE]

for $t$ in some neighborhood of [math], and differentiation yields

[TABLE]

The remaining task is to simplify this expression.

Theorem 1**.**

1. (a)

For every $v\in W$ there exists a unique $R(v)\in\mathbb{R}^{s\times n}$ such that

[TABLE] 2. (b)

The reduced system, in parameterized version (8), is given by

[TABLE] 3. (c)

The matrix $R(v)$ is uniquely determined by the conditions

[TABLE]

and therefore can be obtained from the matrix equation

[TABLE]

with $\left(D\Phi(v)\,|\,P(\Phi(v))\right)$ invertible. In particular, $v\mapsto R(v)$ is smooth. 4. (d)

For every $x\in Z$ let $L(x)\in\mathbb{R}^{s\times 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(\Phi(v))$. Then

[TABLE]

and the reduced system, in parameterized form, is given by

[TABLE]

Proof.

For every $v\in W$ one has $h^{(0)}(\Phi(v))=0$, and by differentiation

[TABLE]

Thus the image of $D\Phi(v)$ is contained in the kernel of $Dh^{(0)}(\Phi(v))$, and these two vector spaces have dimension $s$, hence they are equal. In turn, for $x\in Z$ the kernel of $Dh^{(0)}(x)$ is by construction equal to the image of $Q(x)$. Thus, for every $v$ the matrices $Q(\Phi(v))$ and $D\Phi(v)$ have the same column space, and the latter has full rank. Therefore every column of $Q(\Phi(v))$ is a unique linear combination of the columns of $D\Phi(v)$. Rewritten in matrix language, this is the assertion of part (a). Part (b) is now obvious from equation (8) and injectivity of $D\Phi(v)$.

The first condition given in part (c) is a consequence of part (a), the identity $Q(x)\cdot P(x)=0$ for all $x\in Z$ (which is readily verified from the defining equation (5)), and the fact that $D\Phi(v)$ is an injective linear map. The second condition follows from the fact that $D\Phi(v)\cdot R(v)=Q(\Phi(v))$ is a projection of rank $s$, by using Lemma 1 in the Appendix. Invertibility of the matrix $\left(D\Phi(v)\,|\,P(\Phi(v))\right)$ follows from the direct kernel–image decomposition with respect to $Dh^{(0)}(\Phi(v))$, since the columns of $D\Phi(v)$ span the kernel and the columns of $P(\Phi(v))$ span the image (see [11] for more details).

To prove (d), first recall from Remark 1 that $Dh^{(0)}(x)$ and $P(x)$ have the same column space, therefore $L(x)P(x)=0$ on $Z$. This and $R(v)P(\Phi(v))=0$ from part (c) imply that $R(v)=\Lambda(v)L(\Phi(v))$ for all $v\in W$, with $\Lambda(v)\in\mathbb{R}^{s\times s}$ uniquely determined. Using now the second condition in part (c), we have

[TABLE]

hence $L(\Phi(v))D\Phi(v)$ is invertible and

[TABLE]

which leads to the asserted expression. ∎

Remark 3**.**

1. (a)

To determine the reduced equation (10), there is no need for explicit knowledge of the projection matrix $Q$, or of the matrix $P$ from the decomposition. However, the column space of $Dh^{(0)}(x)$, $x\in Z$, is a crucial ingredient. 2. (b)

On the other hand, knowledge of $P$ and $\mu$ seems indispensable for the computation of $A(x)=D\mu(x)P(x)$, and of $A(\Phi(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**.**

Assume that

[TABLE]

and partition

[TABLE]

Then

[TABLE]

with

[TABLE]

In these expressions the argument of $D\Phi_{1}$ and $D\Phi_{2}$ is $v$ and the argument of $P_{1}$ and $P_{2}$ is $\Phi(v)$.

In the special case when $\Phi_{1}(v)=v$ we get

[TABLE]

Proof.

With $R_{i}$ given as above one verifies

[TABLE]

by direct computation. Rewriting, one obtains

[TABLE]

and this is the defining property of $R$ in Theorem 1. ∎

Remark 4**.**

1. (a)

Up to a relabeling of variables in $\mathbb{R}^{n}$, a partitioning for $\Phi(v)$ as required in Corollary 1 always exists locally, due to the rank condition on the derivative. 2. (b)

The special case $\Phi_{1}(v)=v$ occurs when the critical manifold is the graph of some function. For this case, reduction formulas were derived earlier by Fenichel [9], Lemma 5.4, and Stiefenhofer [21], Equation (2.13). 3. (c)

In the yet more special case that $\Phi_{1}(v)=v$ and $\Phi_{2}(v)=0$ the procedure yields the familiar quasi-steady state reduction. Indeed, in this case one has $\mu(x)=x_{2}$ for $x=(x_{1},\,x_{2})^{\rm tr}$ and $x_{1}\in\mathbb{R}^{s}$, thus $R_{1}=I_{s}$, $R_{2}=-P_{1}P_{2}^{-1}$ and the reduced equation is

[TABLE]

Ignoring higher order terms in $\varepsilon$ (which are irrelevant for Tikhonov-Fenichel reduction) and renaming variables, one obtains the same system by setting the second part of

[TABLE]

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 [12], Proposition 5.)

Finally, with a view on chemical reaction networks, we address conservation laws.

Proposition 2**.**

Let the smooth real-valued function $\psi$ be defined on some open subset of $U$ which has nonempty intersection with $\Phi(W)$, and assume that $\psi$ is a first integral of system (1) for every $\varepsilon$. Then $\widetilde{\psi}:=\psi\circ\Phi$ is constant or a first integral of the parameterized reduced system (9).

Proof.

By [16], Proposition 8 the restriction of $\psi$ to the critical manifold $Z$ is also a first integral of the reduced system (4), thus $D\psi(x)Q(x)h^{(1)}(x)=0$ for all $x\in Z$. Therefore

[TABLE]

which is the characterizing property for first integrals of system (9). ∎

2.3 Illustrative examples

The following small examples have the primary function to illustrate the arguments and reduction procedures from the previous subsection.

We consider a (hypothetical) slow-fast system, with fast reaction

[TABLE]

and slow reaction

[TABLE]

with associated differential equation (according to the procedure from Subsection 3.1 below)

[TABLE]

The critical manifold $Z$ is determined by $Kx_{1}x_{2}=x_{3}$, with $K=k_{1}/k_{-1}$, and we have $P=(1,1,-1)^{\rm tr}$, $\mu(x)=(-k_{1}x_{1}x_{2}+k_{-1}x_{3})$. A parameterization of $Z$ is given by

[TABLE]

hence to determine $R(v)$ via Theorem 1(c) we have to solve

[TABLE]

By straightforward calculations one obtains

[TABLE]

With

[TABLE]

the reduced system becomes

[TABLE]

This is a case where the critical manifold is the graph of the rational function $(x_{1},x_{2})\mapsto Kx_{1}x_{2}$; thus Corollary 1 would also be applicable. Moreover we get

[TABLE]

hence linear stability of the critical manifold follows by Remark 3(b). 2. 2.

As a non-hypothetical variant we discuss the system with the same fast reaction as in part 1, but with slow reaction

[TABLE]

and associated differential equation

[TABLE]

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

[TABLE]

The reduced system becomes

[TABLE]

and we note a further built-in reduction: The differential equation for the reaction network admits the first integral $\psi=x_{1}+x_{3}$ from stoichiometry, hence by Proposition 2 the reduced equation inherits the first integral $\widetilde{\psi}=v_{1}+Kv_{1}v_{2}$. Thus one ends up with a one dimensional reduced equation, as it should be. 3. 3.

For contrast, consider the hypothetical slow-fast system with fast reaction

[TABLE]

and the same slow reaction as in part 1. The differential equation now becomes

[TABLE]

and the critical manifold $Z$ is given by $Kx_{1}^{2}x_{2}^{2}=x_{3}^{3}$, with $K=k_{1}/k_{-1}$, and $P=(2,2,-3)^{\rm tr}$. A parameterization of $Z$ is given by

[TABLE]

It is obvious that

[TABLE]

is of rank two and satisfies $L\cdot P=0$. This yields, by Theorem 1(d),

[TABLE]

Finally the reduced system is given by

[TABLE]

One can further reduce the dimension to one by utilizing the first integral $\psi=4x_{1}+5x_{2}+6x_{3}$ from stoichiometry.

In this example we could have chosen a different parameterization

[TABLE]

for $Z$, which directly represents $Z$ as the graph of a function, but it may be more convenient to work with a reduced system that has rational right hand side. 4. 4.

Finally we sketch an example that is motivated by mechanics, to illustrate that the range of applications does not only include reaction networks. (See Arnold and Anosov [1] for background, and also [24].) Specifically we look at a pair of coupled nonlinear oscillators

[TABLE]

with irrational $\omega>0$, thus we are in a non-resonant scenario. Computing a normal form up to degree three and reduction by invariants $y_{1}^{2}+y_{2}^{2}$, $y_{3}^{2}+y_{4}^{2}$ yields a two-dimensional system, which generically allows to decide about stability. But here we look at a degenerate case, with reduced equation

[TABLE]

and parameters $a<0$, $b>0$ and $c<0$. The choice of signs ensures that the line of stationary points given by $\mu:=ax_{1}+bx_{2}=0$ 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 [math] as $t\to\infty$.

We consider (11) as fast part $h^{(0)}$ of a singularly perturbed system, thus we have

[TABLE]

A straightforward calculation shows that $A(\Phi(v))=ab(1-c)v<0$ whenever $v>0$, so the line $\mu=0$ is attracting in the positive quadrant. Moreover one may choose $L(x)=(cx_{2},\,-x_{1})$ and thus obtains the reduced equation

[TABLE]

for arbitrary small perturbation $h^{(1)}$. The choice

[TABLE]

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.

3 Applications to reaction networks

While the range of applications of Theorem 1 is not restricted to chemical reaction networks, it is natural to discuss these in greater detail: The consideration of slow and fast reactions leads to critical manifolds that consist of stationary points of a subnetwork, and for some relevant and familiar classes of networks explicit parameterizations of the variety of stationary points exist. We first recall some general facts about reaction networks and then discuss two special classes. The results will be illustrated by examples.

3.1 Reaction networks

We briefly recall here the mathematical description of reaction networks according to Feinberg [5], Horn and Jackson [14] (see also the recent monograph [6] by Feinberg), 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},\dots,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=\sum_{i=1}^{n}a_{i}X_{i}$ and is identified with the vector $y=(a_{1},\dots,a_{n})\in\mathbb{R}^{n}$. We let $x_{i}$ denote the concentration of $X_{i}$ and $x=(x_{1},\dots,x_{n})$. For each reaction (denoted by $y\rightarrow y^{\prime}$) we assume given a rate function $w_{y\rightarrow y^{\prime}}(x)\in\mathbb{R}_{\geq 0}$ for $x\in\mathbb{R}^{n}_{\geq 0}$. This leads to a system of differential equations describing the evolution of the concentrations in time:

[TABLE]

Note that $y^{\prime}-y\in\mathbb{R}^{n}$ encodes the net production of each species by the occurrence of the reaction $y\rightarrow y^{\prime}$. The vector subspace spanned by all the $y^{\prime}-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^{\prime}-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

[TABLE]

A frequent choice of rate function is the one from mass action kinetics, with

[TABLE]

with $k_{y\rightarrow y^{\prime}}>0$ called reaction rate constants and using the convention $0^{0}=1$. For the following we recall some definitions:

Definition 1**.**

Let $x,y\in\mathbb{R}^{n}$ and $M\in\mathbb{R}^{n\times m}$, with columns $M_{1},\ldots,M_{m}\in\mathbb{R}^{n}$.

1. (a)

For $x\in\mathbb{R}_{>0}^{n}$ we define

[TABLE]

noting that the definition may be extended to all $x\in\mathbb{R}^{n}$ when all $y_{i}$ are nonnegative integers. 2. (b)

For $x\in\mathbb{R}_{>0}^{n}$ we define

[TABLE]

noting that the definition may be extended to all $x\in\mathbb{R}^{n}$ when all entries of $M$ are nonnegative integers. 3. (c)

The Hadamard product $x\circ y$ is defined as the componentwise multiplication of the two vectors $x,y$, i.e.

[TABLE]

In view of the last definition we can rewrite the reaction network for mass action kinetics in the form

[TABLE]

where $K\in\mathbb{R}_{>0}^{m}$ ($m$ is the number of reactions) is a vector containing the reaction rate constants and $Y\in\mathbb{R}^{n\times 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 [4]. 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 [15] or Feinberg [3].

We turn now to a scenario with prescribed slow and fast reactions; see also the discussions in Heinrich and Schauer [13], Lee and Othmer [17]. 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_{\mathsf{f}}$ and rate vector $w_{\mathsf{f}}(x)$, such that, in the notation of Section 2,

[TABLE]

Analogously, we have

[TABLE]

for the slow subsystem. Keeping the notation from Section 2, we let $Z$ be the zero set of $h^{(0)}$ (possibly restricted to a neighborhood $\widetilde{U}$ of some point), and let $r$ denote the rank of $Dh^{(0)}(x)$, $x\in Z$. Then clearly $\operatorname*{rank}N_{\mathsf{f}}\geq r$, but the inequality 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 when equality holds:

Blanket hypothesis. We impose on system (16) the conditions of Proposition 1(a) and the additional condition that $\operatorname*{rank}N_{\mathsf{f}}={\rm rank}\,Dh^{(0)}(x)=r$, $x\in Z$.

Due to nonnegativity of concentrations, the points of $Z$ will be in $\mathbb{R}^{n}_{\geq 0}$. In some instances we will require that the neighborhood $\widetilde{U}$ in Proposition 1 is even a subset of $\mathbb{R}^{n}_{>0}$, and likewise we will occasionally require that the domain $W$ of the parameterization $\Phi$ is a subset of $\mathbb{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

[TABLE]

with $\operatorname*{rank}D\Phi(v)=s$ on $W$.

Proposition 3**.**

Let system (16) be given, with a parameterization $\Phi$ of the critical manifold as in (7), and assume the blanket hypothesis holds on $\Phi(W)$. Let $L_{\mathsf{f}}\in\mathbb{R}^{s\times n}$ be a matrix whose rows form a basis of the left-kernel of $N_{\mathsf{f}}$. Then the matrix $R(v)$ in Theorem 1 is given as

[TABLE]

and the reduced system is

[TABLE]

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_{\mathsf{f}}$.

Proof.

This follows from Theorem 1(d), since $L_{\mathsf{f}}Dh^{(0)}(x)=L_{\mathsf{f}}N_{\mathsf{f}}Dw_{\mathsf{f}}(x)=0$ on $Z$. ∎

Note that by Remark 2, if $\operatorname*{rank}N=\operatorname*{rank}N_{\mathsf{f}}$, then the reduced system is given by $v^{\prime}=0.$

3.2 Canonical parameterizations for some classes

To find a function $\Phi(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_{\mathsf{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.

3.2.1 Non-interacting sets and rational parameterizations

As was shown in [7], 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 $\mathcal{Y}=\{X_{i_{1}},\dots,X_{i_{r}}\}$, with the following assumptions.

1. (i)

For every fast reaction $y\rightarrow y^{\prime}$, both the sum of the coefficients of the species in $\mathcal{Y}$ in $y$ and the sum of the coefficients in $y^{\prime}$ is at most one. This means that no pair of species in $\mathcal{Y}$ appear together at one side of a reaction, and further no species appears with coefficient greater than $1$. 2. (ii)

The rank of the submatrix of $N_{\mathsf{f}}$ given by the rows $i_{1},\dots,i_{r}$ is equal to $r$. 3. (iii)

Consider the network induced by the fast subnetwork by setting all species not in $\mathcal{Y}$ to zero. For each species $X_{i_{j}}$ in $\mathcal{Y}$, there is a directed path from $X_{i_{j}}$ to [math] in this induced network.

In the nomenclature of [7], assumption (i) means that $\mathcal{Y}$ is non-interacting, (ii) means that no conservation law has support in $\mathcal{Y}$, and (iii) means that there exists a spanning tree rooted at $*$ in the appropriate digraph (see [7], Section 8 for details).

Let $X_{\ell_{1}},\dots,X_{\ell_{s}}$ be the species not in $\mathcal{Y}$. If $\mathcal{Y}$ satisfies (i), (ii) and (iii), then the components $i_{1},\dots,i_{r}$ of $h^{(0)}(x)$ form a linear system in $x_{i_{1}},\dots,x_{i_{r}}$ that has a unique solution in terms of $x_{\ell_{1}},\dots,x_{\ell_{s}}$. Furthermore, the solution is a rational function in $x_{\ell_{1}},\dots,x_{\ell_{s}}$ and in the reaction rate constants $k_{y\rightarrow y^{\prime}}>0$, with all coefficients 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_{\ell_{i}}$, $i=1,\dots,s$. Further, clearly $\operatorname*{rank}D\Phi(v)=s$. In [8] some conditions are stated which guarantee that the assumptions in Proposition 1(a) are satisfied.

3.2.2 Monomial parameterizations and deficiency zero networks

We next consider another common scenario occurring, for instance, for so-called complex balanced steady states (see Feinberg [5], Horn and Jackson [14]) and networks with toric steady states (see Pérez Millán et al. [20], Müller et al. [18]). In this scenario the zero set $Z$ of $h^{(0)}$ in $\mathbb{R}^{n}_{>0}$ agrees with the solution set of a collection of binomial equations

[TABLE]

where $u_{\ell},{c}_{\ell}\in\mathbb{R}^{n}$ and $a_{\ell}(k),b_{\ell}(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

[TABLE]

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\in\mathbb{R}^{n\times{q}}$ the matrix whose columns are $u_{\ell}-{c}_{\ell}$. If $x$ is a solution to (19) then

[TABLE]

It is a classical result (see for example Lemma 3.7 in Müller et al. [18]) that the solution set to this equation, and hence $Z$, can be parameterized in the form

[TABLE]

Here $d=\dim\ker M^{\rm tr}$, and $b_{1},\dots,b_{n}$ are the columns of a matrix $B\in\mathbb{R}^{d\times n}$ with row span equal to $\ker M^{\rm tr}$ and $\ker B^{\rm tr}=\{0\}$ (thus $d={\rm rank}\,B$). With an easy computation one verifies the well-known identity

[TABLE]

where, for a vector $\alpha$, $\operatorname*{diag}(\alpha)$ denotes the diagonal matrix with the entries of $\alpha$ in the diagonal, and $1/v$ is defined component-wise. By this identity, the rank of $D\Phi(v)$ equals $d$, the rank of $B$, and therefore we are in the setting of Subsection 2.2 provided that $d=s$. With this in place, by Proposition 3 the matrix $R(v)$ in Theorem 1 becomes

[TABLE]

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

[TABLE]

As shown in Feinberg [5] and Horn [15], a necessary condition for complex balanced 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\rightarrow y^{\prime}}$ 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

[TABLE]

for every pair of nodes $y_{i},y_{j}$ in the same connected component of the reaction network, and with $K_{ij}$ and $K_{ji}$ positive polynomials in the parameters $k_{y\rightarrow y^{\prime}}$ for the reactions in the same connected component [2]. This implies that the column span of the matrix $M$ above agrees with the column span of $N_{\mathsf{f}}$, and therefore $\ker M^{\rm tr}$ has rank $s$. As a suitable matrix $B$ one can choose $L_{\mathsf{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:

1. (a)

With the notation of Theorem 1 and Proposition 3, there is a parameterization (20) of the critical manifold with $B=L_{\mathsf{f}}$, parameter space $\mathbb{R}_{>0}^{s}$ and

[TABLE] 2. (b)

The reduced system is given by:

[TABLE]

Proof.

Part (a) is clear, while part (b) follows immediately with part (a) and

[TABLE]

∎

3.2.3 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 [15], which are of primary interest in our setting, since we only require positivity of reaction constants in our considerations. By Horn [15], 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 [5] 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 non-zero eigenvalues of the Jacobian have negative real part, hence all the conditions of Proposition 1(a) hold.

3.2.4 Examples

Example 1**.**

We consider the following reaction network with mass action kinetics, where the numbers $k_{y\rightarrow y^{\prime}}$ are written as labels of the reactions:

[TABLE]

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},\dots,k_{6}$, such that the fast subnetwork is

[TABLE]

and the slow reactions have labels $k_{7},k_{8},k_{9}$. With this choice of slow-fast reactions, we have

[TABLE]

and

[TABLE]

The fast reaction network in (26) has $6$ nodes and three connected components, and the rank of $N_{\mathsf{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

[TABLE]

and we easily verify that

[TABLE]

is a positive steady state of the fast subnetwork. With Proposition 4(b) the reduced system can be computed. We choose

[TABLE]

and obtain the following parameterization of $Z$:

[TABLE]

Using equation (24), the reduced system is found:

[TABLE]

where $\xi$ is given by

[TABLE]

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 $\Phi(v)$ is not unique. For example, choosing another starting steady state

[TABLE]

we obtain the parameterization

[TABLE]

and the reduced system

[TABLE]

with

[TABLE]

The same parameterization is obtained by eliminating $x_{2},x_{4},x_{5}$ after realizing that the set of species $\{X_{2},X_{4},X_{5}\}$ satisfies (i), (ii) and (iii) in Subsection 3.2.1.

Example 2**.**

If we remove the reactions with label $k_{8},k_{9}$ from the network (25), then the stoichiometric matrices of both the fast subnetwork and the full network have rank $3$. Hence, if the reduction with a parameterized critical manifold is possible, the reduced system is $v^{\prime}=0$.

Example 3**.**

We analyse the reaction network

[TABLE]

This system can be interpreted as a dual phosphorylation cycle with $X_{1}$ the kinase catalysing the phosphorylation of a substrate $S$ with two phosphorylation sites. Then $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}\leftrightarrow X_{5}$, namely those with label $k_{1},\dots,k_{6}$. Hence the reactions with label $k_{7},k_{8}$ are slow. With this choice, we have

[TABLE]

and

[TABLE]

The fast network has deficiency $1$, since the rank of $N_{\mathsf{f}}$ is $3$, and the network has 6 nodes and two connected components. Thus, the steady states are not complex balanced for all $k$. Instead, we observe that the set $\{X_{3},X_{4},X_{5}\}$ 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 to zero is

\textstyle{0\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{X_{3}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{X_{5}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{X_{4}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}

and clearly there is directed path to [math] from every species. This implies that $x_{3},x_{4},x_{5}$ can be solved from the system $h^{(0)}(x)_{3,4,5}=0$ to obtain the following parameterization (where $v_{1}=x_{1},v_{2}=x_{2},v_{3}=x_{6}$)

[TABLE]

which has an entry that is not monomial. Next we compute the matrix $R(v)$ and the reduced system using Proposition 3. We have

[TABLE]

and using R(v_{1},v_{2},v_{3})=\big{(}L_{\mathsf{f}}\,D\Phi(v)\big{)}^{-1}L_{\mathsf{f}} the reduced system is given by

[TABLE]

where

[TABLE]

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\mu(x)\cdot P(x)$ with

[TABLE]

Using the Routh-Hurwitz conditions (see Gantmacher [10], Ch. V, §6) for its characteristic polynomial $\chi_{A}(\lambda)=\lambda^{3}+\sigma_{1}\lambda^{2}+\sigma_{2}\lambda+\sigma_{3}$ of degree $3$, we obtain:

[TABLE]

In order to verify that all eigenvalues have negative real part we use the Hurwitz conditions for polynomials of degree three:

[TABLE]

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$.

Appendix

For convenient reference, we record a basically known result.

Lemma 1**.**

Let $0<s<n$ and $A\in\mathbb{R}^{n\times s}$, $B\in\mathbb{R}^{s\times n}$ such that ${\rm rank}\,AB=s\text{ and }(AB)^{2}=AB.$ Then $BA=I_{s}$.

Proof.

Recall that $\mathbb{R}^{n}$ is the direct sum of the eigenspaces of $AB$ for eigenvalues $1$ and [math]. Let $z_{1},\ldots,z_{s}$ be a basis of the former. Then $Bz_{1},\ldots,Bz_{s}$ are linearly independent due to

[TABLE]

hence they form a basis of $\mathbb{R}^{s}$. Finally, $ABz_{i}=z_{i}$ implies $(BA)Bz_{i}=Bz_{i}$ for $1\leq i\leq s$. ∎

Acknowledgements. The work of NK and SW has been supported by the bilateral project ANR-17-CE40-0036 and DFG-391322026 SYMBIONT. EF has been supported by the Independent Research Fund of Denmark. SW thanks the Department of Mathematical Sciences of the University of Copenhagen for its hospitality during a research visit when essential parts of the present manuscript were devised. Likewise, EF thanks the hospitality of RWTH Aachen where this project was initiated.