Power-engine-load form for
dynamic absolute concentration robustness

Badal Joshi ¹¹1Department of Mathematics, California State University San Marcos (bjoshi

@

csusm.edu) and Gheorghe Craciun ²²2Departments of Mathematics and Biomolecular Chemistry, University of Wisconsin-Madison (craciun

@

wisc.edu)

Abstract

In a reaction network, the concentration of a species with the property of dynamic absolute concentration robustness (dynamic ACR) converges to the same value independent of the overall initial values. This property endows a biochemical network with output robustness and therefore is essential for its functioning in a highly variable environment. It is important to identify structure of the dynamical system as well as constraints required for dynamic ACR. We propose a power-engine-load form of dynamic ACR and obtain results regarding convergence to the ACR value based on this form.
Keywords: biochemical reaction network, absolute concentration robustness, power-engine-load form, dynamical systems, robust network output, chemostat.

AMS subject codes: 34C20, 37N25, 37N35, 92C42

1 Introduction

Dynamic absolute concentration robustness (dynamic ACR), introduced in [1], is a property of dynamical systems wherein one variable converges to a unique value independent of the initial values. This variable is the dynamic ACR variable and the unique value is its ACR value. Dynamic ACR is significant for applications to biochemistry. Biochemical systems need to perform robustly in a wide variety of conditions. Dynamic ACR provides a mechanism for such robustness. In signal response circuits, an essential design feature is that the response depends on the signal but not on the internal state of the signaling circuit. Different initial states (which model the internal state of the signaling circuit when the signal arrives) generically result in different final states, especially when some conservation laws hold. However, if one of the variables is invariant across all final states, such as is the case in dynamic ACR, then this variable can be taken to be the invariant signal response. The dependence of the signal strength is encoded in the reaction rate constants, and therefore signal response will depend on the signal strength via the rate constants. However, the signal response will not depend on the initial values, i.e. the internal state of the signaling circuit.

Formally, dynamic ACR is the property that all trajectories (with some minor restrictions on allowed initial values) converge to the hyperplane $\{x_{i}=a_{i}^{*}\}$ . It is desirable to identify sufficient conditions, either on the underlying reaction network or the dynamical system, that generate dynamic ACR. Necessary and sufficient conditions for dynamic ACR in complex balanced systems have been obtained in [1]. Moreover, a classification of small networks with dynamic ACR appears in [2].

Shinar and Feinberg [3] defined static absolute concentration robustness (static ACR) as the property that all positive steady states are located in the hyperplane $\{x_{i}=a_{i}^{*}\}$ . Furthermore, they gave sufficient network conditions for the reaction network to have static ACR in some species. Additionally, they gave many biochemically realistic networks and showed that these networks satisfy their network conditions and therefore have static ACR. Many of the biochemical networks appearing in their paper may have the property of dynamic ACR as well, although this remains an open question. It is difficult to characterize dynamic ACR because it requires understanding the limiting behavior for arbitrary initial values. We propose a power-engine-load form (see (3.1)) for the differential equation satisfied by the supposed dynamic ACR variable, which may help establish dynamic ACR in some situations as we show in this paper. In Theorem 3.8, we establish sufficient conditions for convergence to the ACR value, based on properties of the terms appearing in the power-engine-load form equation.

Dynamic ACR has been experimentally observed in bacterial two-component signaling systems such as the EnvZ-OmpR system and the IDHKP-IDH system [4, 5, 6, 7, 8, 9]. Absolute concentration robustness (both static and dynamic) is related to the concept called robust perfect adaptation [10, 11, 12] studied from a control theory perspective. For a biochemical perspective on perfect adaptation, see [13, 14]. Structural requirements for robust perfect adaptation in biomolecular networks are studied in [15]. While our goal in this paper is not to explore the connections between the two fields of dynamical systems and control theory, we mention the above papers here as a help to any readers interested in understanding the connection further.

This article is organized as follows. Section 2 gives the background on reaction networks, mass action kinetics, and dynamic absolute concentration robustness (ACR). Section 3 introduces a power-engine-load form for dynamic ACR, examples, and convergence results. The main result of this paper is Theorem 3.8. Section 4 studies the dynamics of several reaction networks where we apply power-engine-load form and illustrate the use of Theorem 3.8. In many cases, the reaction network is coupled with inflows. While each network requires special analysis, we ultimately apply Theorem 3.8 to show dynamic ACR for the network under study.

2 Reaction networks

An example of a reaction is $A+B\to 2C$ , which is a schematic representation of the process where a molecule of species $A$ and a molecule of species $B$ react with one another and result in two molecules of a species $C$ . The abstract linear combination of species $A+B$ that appears on the left of the reaction arrow is called the reactant complex while $2C$ is called the product complex of the reaction $A+B\to 2C$ . We assume that for any reaction the product complex is different from the reactant complex. A reaction network is a nonempty, finite set of species and a nonempty, finite set of reactions such that every species appears in at least one complex. We will use mass action kinetics for the rate of reactions.

2.1 Mass action systems

In mass action kinetics, each reaction occurs at a rate proportional to the product of the concentrations of species appearing in the reactant complex. We conventionally use lower case letters $a,b,c$ to denote the species concentrations of the corresponding species $A,B,C$ , respectively. Under mass action kinetics, the reaction $A+B\to 2C$ occurs at rate $kab$ where $k$ is the reaction rate constant, conventionally placed near the reaction arrow, as follows: $A+B\xrightarrow{k}2C$ . Consider the reaction network $\{A+B\xrightarrow{k_{1}}2C,~{}2C\xrightarrow{k_{2}}2A\}$ . Application of mass action kinetics results in the following dynamical system, called mass action system, that describes the time-evolution of the species concentrations.

	$\displaystyle\dot{a}$	$\displaystyle=-k_{1}ab+2k_{2}c^{2}$
	$\displaystyle\dot{b}$	$\displaystyle=-k_{1}ab$
	$\displaystyle\dot{c}$	$\displaystyle=2k_{1}ab-2k_{2}c^{2}$

For further details on reaction networks and mass action systems, see for instance [16, 17]. We use standard notation and terminology of dynamical systems, such as steady states, stability, basin of attraction etc., see for instance [18, 19].

2.2 Dynamic absolute concentration robustness

Shinar and Feinberg defined (static) absolute concentration robustness (static ACR) as the property that all positive steady states are contained in a fixed hyperplane parallel to a coordinate hyperplane [3]. This means that for every positive steady state, one of the coordinates is invariant. Several network/structural conditions for static ACR have been identified [3, 20, 1, 2]. Static ACR is designed to model the property of a robust signal response despite variability in the internal state of the signaling circuit. However, the property accurately models output robustness only if there is convergence to a positive steady state for every initial value, and indeed, such convergence is not required in the definition of static ACR.

We defined dynamic absolute concentration robustness (dynamic ACR) as the property that all initial values lead to convergence to a fixed hyperplane parallel to a coordinate hyperplane [1] – a requirement somewhat more general than converging to a steady state in the hyperplane. We believe that this property better captures the idea of output robustness, when compared to static ACR, since it models long-term dynamics. Network conditions for dynamic ACR in small reaction networks were found in [2]. The condition found there is geometric in nature – for dynamic ACR in networks with two reactions and two species, the reactant polytope (line segment joining the two reactant complexes in their geometric embedding) must be parallel to a coordinate axis. We seek an extension of this condition in arbitrary networks. If the differential equation for the candidate ACR species has a certain special form (power-engine-load form), then it may have dynamic ACR provided that the power and load satisfy some additional conditions.

Dynamic ACR was defined for autonomous systems in [1]. We generalize to include non-autonomous systems. Throughout the paper, we assume that $\mathcal{D}$ is a dynamical system defined by $\dot{x}=F(x,t)$ with $x\in\mathbb{R}^{n}_{\geq 0}$ for which $\mathbb{R}^{n}_{\geq 0}$ is forward invariant.

Definition 2.1.

The kinetic subspace of $\mathcal{D}$ is defined to be the linear span of the image of $F$ , denoted by ${\rm span}(\imaginary(F))$ . Two points $x,y\in\mathbb{R}^{n}_{\geq 0}$ are compatible if $y-x\in{\rm span}(\imaginary(F))$ . The sets $S,S^{\prime}\subseteq\mathbb{R}^{n}_{\geq 0}$ are compatible if there are $x\in S$ and $x^{\prime}\in S^{\prime}$ such that $x$ and $x^{\prime}$ are compatible. A compatibility class $S$ is a nonempty subset of $\mathbb{R}^{n}_{\geq 0}$ such that $x,y\in S$ if and only if $y-x\in{\rm span}(\imaginary(F))$ .

Definition 2.2.

$\mathcal{D}$ is a dynamic ACR system if there is an $i\in\{1,\ldots,n\}$ with $F_{i}\not\equiv 0$ and a positive $a_{i}^{*}\in\mathbb{R}_{>0}$ such that for any $(t_{0},x(t_{0}))\in\mathbb{R}\times\mathbb{R}^{n}_{>0}$ with $x(t_{0})$ is compatible with $\{x\in\mathbb{R}^{n}_{>0}~{}|~{}x_{i}=a_{i}^{*}\}$ , a unique solution to $\dot{x}=F(x,t)$ exists up to some maximal $T_{0}(t_{0},x(t_{0}))\in(t_{0},\infty]$ , and $x_{i}(t)\xrightarrow{t\to T_{0}}a_{i}^{*}$ . Any such $x_{i}$ and $a_{i}^{*}$ is a dynamic ACR variable and its dynamic ACR value, respectively.

If the dynamical system $\dot{x}=F(x,t)$ does not have the possibility of a finite-time blow-up, then $T_{0}(t_{0},x(t_{0}))=\infty$ for any $(t_{0},x(t_{0}))\in\mathbb{R}\times\mathbb{R}^{n}_{>0}$ . In this paper, we focus on dynamic ACR in systems where a unique solution is assumed to exist for all positive time.

3 Power-engine-load form of dynamic ACR

In [2], we discuss the network conditions for static and dynamic ACR in reaction networks with two species and two reactions. One of the conditions is that the reactant polytope (i.e. in the Euclidean embedding of the reaction network, the line segment joining the two reactant complexes) is parallel to one of the coordinate axes. A natural generalization of this geometric condition is the one given below, called power-engine-load form for dynamic ACR.

Suppose that $x\in\mathbb{R}^{n}$ and for some $1\leq i\leq n$ , $x_{i}$ satisfies

\displaystyle\derivative{x_{i}}{t}~{}=~{}{\color[rgb]{0.5,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.5,0.5,0}\pgfsys@color@cmyk@stroke{0}{0}{1}{.5}\pgfsys@color@cmyk@fill{0}{0}{1}{.5}\underbrace{f(x(t),t)}_{\mbox{power}}}~{}\cdot~{}{\color[rgb]{1,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,1}\pgfsys@color@cmyk@stroke{0}{1}{0}{0}\pgfsys@color@cmyk@fill{0}{1}{0}{0}\underbrace{(x_{i}^{*}-x_{i}(t))}_{\mbox{engine}}}~{}+~{}{\color[rgb]{0,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,1}\pgfsys@color@cmyk@stroke{1}{0}{0}{0}\pgfsys@color@cmyk@fill{1}{0}{0}{0}\underbrace{g(x(t),t)}_{\mbox{load}}}

(3.1)

then under some reasonable conditions on “power” $f(x(t),t)$ and “load” $g(x(t),t)$ , the variable $x_{i}$ will have dynamic ACR with the value $x_{i}^{*}$ . An example of sufficient conditions is $g\equiv 0$ , $f>0$ , and $\int_{0}^{\infty}f(x(t),t)~{}dt=\infty$ . We state the precise result and prove other convergence results in Section 3.2.

3.1 Examples of power-engine-load form in mass action systems

Example 3.1.

The simplest example of a mass action system with power-engine-load form is the reaction network $0\stackrel{{\scriptstyle[}}{{k}}^{\prime}]{k}{\rightleftarrows}X$ and its associated ODE

\derivative{x}{t}={\color[rgb]{0.5,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.5,0.5,0}\pgfsys@color@cmyk@stroke{0}{0}{1}{.5}\pgfsys@color@cmyk@fill{0}{0}{1}{.5}k^{\prime}}{\color[rgb]{1,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,1}\pgfsys@color@cmyk@stroke{0}{1}{0}{0}\pgfsys@color@cmyk@fill{0}{1}{0}{0}\left(k/k^{\prime}-x\right)}.

It is straightforward to show that $x$ has dynamic ACR with value $k/k^{\prime}$ . ∎

Example 3.2.

Biologically interesting cases require a reactant complex with more than one species (so that a ‘reaction’ can occur) and some positive mass conservation law involving all species. The following well-known network is the simplest that satisfies these requirements of species interaction and mass conservation:

\displaystyle A+B\xrightarrow{k_{1}}2B,\quad B\xrightarrow{k_{2}}A.

(3.2)

The positive mass conservation law $a+b=const$ is apparent from the mass action ODE system:

\displaystyle\dot{a}={\color[rgb]{0.5,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.5,0.5,0}\pgfsys@color@cmyk@stroke{0}{0}{1}{.5}\pgfsys@color@cmyk@fill{0}{0}{1}{.5}k_{1}b}{\color[rgb]{1,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,1}\pgfsys@color@cmyk@stroke{0}{1}{0}{0}\pgfsys@color@cmyk@fill{0}{1}{0}{0}\left(k_{2}/k_{1}-a\right)},\quad\dot{b}=-k_{1}b\left(k_{2}/k_{1}-a\right).

(3.3)

The ODE for $a$ has power-engine-load form and it is known (see for instance [1]) that $a(t)\xrightarrow{t\to\infty}k_{2}/k_{1}$ for any initial value $(a(0),b(0))=(a_{0},b_{0})$ that satisfies $a_{0}+b_{0}\geq k_{2}/k_{1}$ . ∎

Example 3.3.

In bacterial two-component signaling systems, the circuit mechanism for robust signal transduction from the cell environment to its interior involves a bifunctional component, a mechanism that is found in thousands of biological systems [8]. One such system is the E. coli IDHKP-IDH glyoxylate bypass regulation system whose core ACR module is

	$\displaystyle X+E$	$\displaystyle\stackrel{{\scriptstyle[}}{{k}}_{2}]{k_{1}}{\rightleftarrows}C_{1}\xrightarrow{k_{3}}Y+E$
	$\displaystyle Y+C_{1}$	$\displaystyle\stackrel{{\scriptstyle[}}{{k}}_{5}]{k_{4}}{\rightleftarrows}C_{2}\xrightarrow{k_{6}}X+C_{1}.$		(3.4)

It is known that (3.3) has static ACR in $Y$ [3]. The mass action ODE equation for the concentration of $Y$ can be written in power-engine-load form, using the fact that the static ACR value of $y$ is $k_{3}/k_{4}(1+k_{5}/k_{6})$ , as follows.

	$\displaystyle\dot{y}$	$\displaystyle=k_{3}c_{1}-k_{4}c_{1}y+k_{5}c_{2}$
		$\displaystyle=k_{3}c_{1}-k_{4}c_{1}y+k_{5}\frac{k_{3}}{k_{6}}c_{1}+\frac{k_{5}}{k_{6}}\left(k_{6}c_{2}-k_{3}c_{1}\right)$
		$\displaystyle={\color[rgb]{0.5,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.5,0.5,0}\pgfsys@color@cmyk@stroke{0}{0}{1}{.5}\pgfsys@color@cmyk@fill{0}{0}{1}{.5}k_{4}c_{1}}{\color[rgb]{1,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,1}\pgfsys@color@cmyk@stroke{0}{1}{0}{0}\pgfsys@color@cmyk@fill{0}{1}{0}{0}\left(\frac{k_{3}}{k_{4}}\left(1+\frac{k_{5}}{k_{6}}\right)-y\right)}+{\color[rgb]{0,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,1}\pgfsys@color@cmyk@stroke{1}{0}{0}{0}\pgfsys@color@cmyk@fill{1}{0}{0}{0}\frac{k_{5}}{k_{6}}\left(k_{6}c_{2}-k_{3}c_{1}\right)}.$

The load is not identically zero in this case. We will prove in future work that $y$ has dynamic ACR with ACR value $k^{*}=\dfrac{k_{3}}{k_{4}}\left(1+\dfrac{k_{5}}{k_{6}}\right)$ . As shown in Theorem 3.8, a sufficient condition for dynamic ACR is that for any positive initial value $\left(k_{6}c_{2}(t)-k_{3}c_{1}(t)\right)\xrightarrow{t\to\infty}0$ and that $\int_{0}^{\infty}c_{1}(t)~{}dt=\infty$ . Checking that the two conditions are satisfied requires some additional ideas which will be developed in future work. ∎

We now discuss conditions on “power” $f(x(t),t)$ and “load” $g(x(t),t)$ that either ensure or prevent dynamic ACR.

3.2 Convergence results for power-engine-load form dynamic ACR

Theorem 3.4 (Zero load).

Consider the dynamical system $\dot{x}=F(x,t)$ with continuously differentiable $F:\mathbb{R}^{n}_{\geq 0}\times\mathbb{R}_{\geq 0}\to\mathbb{R}^{n}$ and for which $\mathbb{R}^{n}_{\geq 0}$ is forward invariant. Suppose there is an $x_{i}^{*}\in\mathbb{R}_{>0}$ such that $F_{i}(x,t)|_{x\in\mathbb{R}^{n}_{>0}}=0$ if and only if $x_{i}=x_{i}^{*}$ . The following hold for every solution $(x(t))_{t\geq 0}$ of $\dot{x}=F(x,t)$ .

1.

$x_{i}^{*}-x_{i}(t)$ has the same sign as $x_{i}^{*}-x_{i}(0)$ for all $t\geq 0$ .
2.

Suppose $x_{i}(0)\neq x_{i}^{*}$ . Then $\displaystyle\frac{F_{i}(x(t),t)}{x_{i}^{*}-x_{i}(t)}$ has the same sign as $\displaystyle\frac{F_{i}(x(0),0)}{x_{i}^{*}-x_{i}(0)}$ for all $t\geq 0$ .

Suppose $x_{i}(0)\neq x_{i}^{*}$ .

(a)

If $\displaystyle\frac{F_{i}(x(0),0)}{x_{i}^{*}-x_{i}(0)}>0$ then $\absolutevalue{x_{i}(t)-x_{i}^{*}}$ is strictly decreasing on $[0,\infty)$ and

\absolutevalue{x_{i}(\infty)-x_{i}^{*}}=\absolutevalue{x_{i}(0)-x_{i}^{*}}\exp\left(-\int_{0}^{\infty}\frac{F_{i}(x(s),s)}{x_{i}^{*}-x_{i}(s)}ds\right).

(b)

If $\displaystyle\frac{F_{i}(x(0),0)}{x_{i}^{*}-x_{i}(0)}<0$ then $\absolutevalue{x_{i}(t)-x_{i}^{*}}$ is strictly increasing on $[0,\infty)$ .

4.

Suppose $x_{i}(0)\neq x_{i}^{*}$ . Then $x_{i}(t)\xrightarrow{t\to\infty}x_{i}^{*}$ if and only if

$\int_{0}^{\infty}\frac{F_{i}(x(t),t)}{x_{i}^{*}-x_{i}(t)}dt=\infty.$

Proof.

If $x_{i}(0)=x_{i}^{*}$ then $\dot{x}_{i}|_{t=0}=F_{i}(x(0),0)=0$ and so $x_{i}(t)=x_{i}^{*}$ for all $t\geq 0$ . If there is a $t>0$ such that $x_{i}(t)-x_{i}^{*}$ and $x_{i}(0)-x_{i}^{*}$ have different signs then by continuity of $x_{i}(t)$ , there must be a $t^{\prime}\in(0,t)$ such that $x_{i}(t^{\prime})=x_{i}^{*}$ . But then $x_{i}(t)=x_{i}^{*}$ for all $t\in\mathbb{R}$ , which is a contradiction.

If $x_{i}(0)\neq x_{i}^{*}$ then by the previous part $\displaystyle\frac{F_{i}(x(t),t)}{x_{i}^{*}-x_{i}(t)}$ is defined for all time $t\geq 0$ . Since $F_{i}(x,t)\neq 0$ for $x_{i}\neq x_{i}^{*}$ , the second result follows. Since $x_{i}(t)\neq x_{i}^{*}$ for all $t\geq 0$ , we can divide by $x_{i}^{*}-x_{i}$ and integrate to get

\displaystyle\frac{dx_{i}}{x_{i}^{*}-x_{i}}=\frac{F_{i}(x(t),t)}{x_{i}^{*}-x_{i}(t)}dt\implies\frac{x_{i}^{*}-x_{i}(t)}{x_{i}^{*}-x_{i}(0)}=\exp\left(-\int_{0}^{t}\frac{F_{i}(x(s),s)}{x_{i}^{*}-x_{i}(s)}ds\right).

(3.5)

By the previous result, the integrand has the same sign as $\displaystyle\frac{F_{i}(x(0),0)}{x_{i}^{*}-x_{i}(0)}$ for all positive time, and so result 3 follows. Finally, result 4 follows from taking $\lim_{t\to\infty}$ on both sides. ∎

Example 3.5.

Consider the mass action system (3.3), $\dot{a}=k_{1}b\left(k-a\right),~{}\dot{b}=-k_{1}b\left(k-a\right)$ . The variable $a$ has power-engine-load form with zero load and power $k_{1}b>0$ for every $(a,b)\in\mathbb{R}^{2}_{>0}$ . It is clear that $a(0)=k$ if and only if $\dot{a}|_{t\geq 0}=0$ . So assume that $a(0)\neq k$ .

The mass action system with initial value $(a(0),b(0))=(a_{0},b_{0})\in\mathbb{R}^{2}_{\geq 0}$ can be solved explicitly after using the conservation relation $a(t)+b(t)=a_{0}+b_{0}$ , thus reducing to the one-dimensional system $\dot{b}=-k_{1}b(k+b-a_{0}-b_{0})$ . We get

\displaystyle b(t)=\begin{cases}\dfrac{a_{0}+b_{0}-k}{1+\left(\frac{a_{0}-k}{b_{0}}\right)e^{-(a_{0}+b_{0}-k)t}}&\mbox{ if }a_{0}+b_{0}\neq k,\\ \dfrac{k-a_{0}}{1+(k-a_{0})t}&\mbox{ if }a_{0}+b_{0}=k.\end{cases}

and $a(t)=a_{0}+b_{0}-b(t)$ .

For $(a_{0},b_{0})\in\mathbb{R}^{2}_{>0}$ , it is well-known that $a(t)\xrightarrow{t\to\infty}k$ if and only if $a_{0}+b_{0}\geq k$ . We now argue that $\int_{0}^{\infty}b(t)~{}dt=\infty$ if and only if $a_{0}+b_{0}\geq k$ and $a_{0}\neq k$ . Indeed, for any $a_{0}<k$ ,

\int_{0}^{\infty}\frac{k-a_{0}}{1+(k-a_{0})t}~{}dt=\infty.

Moreover

\displaystyle\lim_{t\to\infty}\frac{a_{0}+b_{0}-k}{1+\left(\frac{a_{0}-k}{b_{0}}\right)e^{-(a_{0}+b_{0}-k)t}}=\begin{cases}a_{0}+b_{0}-k&\mbox{ if }a_{0}+b_{0}>k,\\ 0&\mbox{ if }a_{0}+b_{0}<k,\end{cases}

(3.6)

implies that

\int_{0}^{\infty}\frac{a_{0}+b_{0}-k}{1+\left(\frac{a_{0}-k}{b_{0}}\right)e^{-(a_{0}+b_{0}-k)t}}~{}dt\begin{cases}=\infty&\mbox{ if }a_{0}+b_{0}>k,\\ <\infty&\mbox{ if }a_{0}+b_{0}<k,\end{cases}

because in the first case ( $a_{0}+b_{0}>k$ ) the integrand does not go to zero, the integral is clearly divergent and in the latter case ( $a_{0}+b_{0}<k$ ) the integrand is $\approx\exp{(a_{0}+b_{0}-k)t}$ , so the integral is convergent.

This shows that $a\xrightarrow{t\to\infty}k$ if and only if either $a(0)=k$ or $\int_{0}^{\infty}k_{1}b(t)~{}dt=\infty$ . ∎

Corollary 3.6.

Suppose the hypotheses of Theorem 3.4 hold and $x_{i}(0)\neq x_{i}^{*}$ . If

\displaystyle\liminf_{t\to\infty}\frac{tF_{i}(x(t),t)}{x_{i}^{*}-x_{i}(t)}\in(0,\infty],

(3.7)

then $x_{i}(t)\xrightarrow{t\to\infty}x_{i}^{*}$ .

Proof.

By (3.7), there is a $\delta>0$ and a $t_{0}\geq 0$ such that $\dfrac{tF_{i}(x(t),t)}{x_{i}^{*}-x_{i}(t)}>\delta$ for all $t\geq t_{0}$ . Then for any $T\geq t_{0}$ ,

\displaystyle\int_{0}^{\infty}\frac{F_{i}(x(t),t)}{x_{i}^{*}-x_{i}}dt\geq\int_{t_{0}}^{T}\frac{F_{i}(x(t),t)}{x_{i}^{*}-x_{i}}dt=\int_{t_{0}}^{T}\frac{1}{t}~{}\frac{tF_{i}(x(t),t)}{x_{i}^{*}-x_{i}}dt>\delta\int_{t_{0}}^{T}\frac{1}{t}dt\xrightarrow{T\to\infty}\infty,

and so by Theorem 3.4, $x_{i}\to x_{i}^{*}$ . ∎

Example 3.7.

The condition (3.7) is sufficient but not necessary for $x_{i}\to x_{i}^{*}$ . Consider the mass action system of the following reaction network:

	$\displaystyle 2X_{1}\xrightarrow{k_{1}}X_{1},$	$\displaystyle\quad 2X_{2}\xrightarrow{k_{2}}X_{2},$
	$\displaystyle X_{1}+X_{2}\xrightarrow{k_{2}}X_{2},$	$\displaystyle\quad Y+X_{1}\stackrel{{\scriptstyle[}}{{k}}_{4}]{k_{5}}{\rightleftarrows}X_{1},$		(3.8)

where the rate constants are the same for the reactions $2X_{2}\to X_{2}$ and $X_{1}+X_{2}\to X_{2}$ ( $X_{2}$ degrades both $X_{1}$ and $X_{2}$ at the same rate). The mass action ODEs are

$\displaystyle\dot{x}_{2}$	$\displaystyle=-k_{2}x_{2}^{2},$
$\displaystyle\dot{x}_{1}$	$\displaystyle=-k_{1}x_{1}^{2}-k_{2}x_{1}x_{2},$
$\displaystyle\dot{y}$	$\displaystyle=k_{5}x_{1}\left(\frac{k_{4}}{k_{5}}-y\right).$	(3.9)

It is simple to check that given an arbitrary positive initial value $(y(0),x_{1}(0),x_{2}(0))=(b_{0},b_{1},b_{2})\in\mathbb{R}^{3}_{>0}$ , the unique solution satisfies

	$\displaystyle x_{2}(t)$	$\displaystyle=\frac{b_{2}}{1+k_{2}b_{2}t},$
	$\displaystyle x_{1}(t)$	$\displaystyle=\frac{b_{1}b_{2}k_{2}}{(1+b_{2}k_{2}t)(b_{2}k_{2}+b_{1}k_{1}\log(1+b_{2}k_{2}t)}\sim\frac{1}{k_{1}t\log(1+b_{2}k_{2}t)}$

and so $\int_{t_{0}}^{\infty}x_{1}(t)dt=\infty$ for any $t_{0}>0$ which implies that $y\to k_{4}/k_{5}$ by Theorem 3.4. However, $\liminf_{t\to\infty}(tx_{1}(t))=0$ , so Corollary 3.6 does not apply. ∎

When the load $g$ is nonzero, the power $f$ must overpower the load $g$ for convergence of $x_{i}$ to $x_{i}^{*}$ .

Theorem 3.8.

Consider the dynamical system $\dot{x}=F(x,t)$ with $F:\mathbb{R}^{n}_{\geq 0}\times\mathbb{R}_{\geq 0}\to\mathbb{R}^{n}$ and for which $\mathbb{R}^{n}_{\geq 0}$ is forward invariant. Suppose that for some $i\in\{1,\ldots,n\}$ we have

F_{i}(x,t)=f(x,t)\cdot\left(x_{i}^{*}-x_{i}\right)+g(x,t),

with $f(x,t)>0$ in $\mathbb{R}^{n}_{>0}\times\mathbb{R}_{>0}$ , and $g\not\equiv 0$ . Let $(x(t))_{t\geq 0}$ be a solution of $\dot{x}=F(x,t)$ such that $\displaystyle\int_{0}^{\infty}f(x(t),t)dt=\infty$ , and such that the limit $\displaystyle\alpha\coloneqq\lim_{t\to\infty}\frac{g(x(t),t)}{f(x(t),t)}$ exists, with $\alpha>-x_{i}^{*}$ . Then we have $x_{i}(t)\xrightarrow{t\to\infty}x_{i}^{*}+\alpha$ .

Proof.

Note that the equation $\displaystyle\derivative{x_{i}}{t}=F_{i}(x,t)$ can be rewritten as

\displaystyle\derivative{x_{i}}{t}

\displaystyle=f(x,t)\left(x_{i}^{*}+\frac{g(x,t)}{f(x,t)}-x_{i}\right).

(3.10)

Then, if we denote $\widetilde{x}_{i}^{*}=x_{i}^{*}+\alpha$ and $\widetilde{g}=g/f-\alpha$ , we can reduce our problem to showing that if $x_{i}(t)$ satisfies the equation

\displaystyle\derivative{x_{i}}{t}

\displaystyle=f(x,t)\left(\ \widetilde{x}_{i}^{*}-x_{i}+\widetilde{g}(x,t)\ \right),

(3.11)

and $\displaystyle\lim_{t\to\infty}\widetilde{g}(x(t),t)=0$ , then $x_{i}(t)\xrightarrow{t\to\infty}\widetilde{x}_{i}^{*}$ .

For any fixed $\varepsilon\in(0,\widetilde{x}_{i}^{*})$ we will now show that there exists some $T_{0}>0$ such that $x_{i}(t)\in(\widetilde{x}_{i}^{*}-\varepsilon,\widetilde{x}_{i}^{*}+\varepsilon)$ for all $t>T_{0}$ . For this, let us first choose some $T_{1}>0$ such that $\widetilde{g}(x(t),t)\in(-\frac{\varepsilon}{2},\frac{\varepsilon}{2})$ for all $t>T_{1}$ . Note that this implies that the interval $(\widetilde{x}_{i}^{*}-\varepsilon,\widetilde{x}_{i}^{*}+\varepsilon)$ is an invariant set of (3.10) for $t>T_{1}$ . Therefore, if $x_{i}(T_{1})\in(\widetilde{x}_{i}^{*}-\varepsilon,\widetilde{x}_{i}^{*}+\varepsilon)$ , we can choose $T_{0}=T_{1}$ . Assume now that $x_{i}(T_{1})\notin(\widetilde{x}_{i}^{*}-\varepsilon,\widetilde{x}_{i}^{*}+\varepsilon)$ , and for example $x_{i}(T_{1})\geq\widetilde{x}_{i}^{*}+\varepsilon$ (the case where $x_{i}(T_{1})\leq\widetilde{x}_{i}^{*}-\varepsilon$ is analogous).

Let us assume that the inequality $x_{i}(t)\geq\widetilde{x}_{i}^{*}+\varepsilon$ holds for all $t>T_{1}$ ; we will show that this leads to a contradiction. Indeed, for any such $t$ we have

\widetilde{x}_{i}^{*}-x_{i}+\widetilde{g}(x,t)<-\frac{\varepsilon}{2},

which implies that

\displaystyle\derivative{x_{i}}{t}

\displaystyle\leq-\frac{\varepsilon}{2}f(x,t),

(3.12)

for all $t>T_{1}$ ; but this, together with the hypothesis that $\displaystyle\int_{0}^{\infty}f(x(t),t)dt=\infty$ , would imply that $\displaystyle\lim_{t\to\infty}x_{i}(t)=-\infty$ , which contradicts our assumption that for all $t>T_{1}$ we have $x_{i}(t)\geq\widetilde{x}_{i}^{*}+\varepsilon$ .

Therefore we obtain the desired conclusion that $x_{i}(t)$ will enter the interval $(\widetilde{x}_{i}^{*}-\varepsilon,\widetilde{x}_{i}^{*}+\varepsilon)$ in some finite time $T_{0}$ , and will remain inside it for all $t>T_{0}$ . ∎

4 Applications of power-engine-load form

We use the power-engine-load form and Theorem 3.8 to prove dynamic ACR in certain reaction networks taken with inflows. The underlying reaction networks have dynamic ACR when there are no inflows or outflows, i.e. when the system is closed. We show that even when certain inflows are included, the property of dynamic ACR persists. Moreover, the ACR value for the open system (with inflows) remains the same in many instances as the ACR value for the closed/isolated system. In Section 4.4, we also consider an application where both inflows and outflows are included, and again given the right conditions the dynamic ACR property persists. In Section 4.7, we consider a simple enzyme catalysis network with a bifunctional enzyme. We show once again that dynamic ACR persists under many different possible inflow conditions. Each individual network requires special analysis in combination with application of Theorem 3.8.

The reaction networks that we consider here have a conserved quantity when the system is closed. When inflows (but not outflows) are included, the total concentration grows in an unbounded manner. This behavior can be realized in a chemostat (for some finite time) but there are also biologically realistic situations where this occurs. For instance, specific ion channels on the surface of a cell might be open and selectively permeable which results in influx of channel-specific ions from outside the cell volume. Our analysis shows that there are invariant quantities even while influxes are ongoing, that is even when the overall system is not near any long-term equilibrium. One variable, the ACR variable, can consistently converge to a robust value even when the overall trajectory does not converge. Robust convergence of one variable when the overall system is far from equilibrium can even occur when inflows and outflows are both present, as shown in Section 4.4. In all these cases, Theorem 3.8 is used to show convergence.

4.1 Motifs of static and dynamic ACR

In [2], we identified the network motifs with two reactions and two species which have both static ACR and dynamic ACR when treated as isolated or closed systems. Each motif represents an entire infinite set of reaction networks. We consider three such motifs shown in Figure 1 and we select one network for each motif. It is shown in [2] that these three are the only motifs (up to interchange of species labels) which have two reactions, two species, and have both properties of static ACR and dynamic ACR.

Figure 1: Motifs with two reactions and two species that have both static and dynamic ACR in one species

A

In each network motif, the arrows represent reactions while the (dashed) line segment joining the two arrow tails is the so-called reactant polytope which plays an important role in the classification of ACR systems [2] – the reactant polytope is parallel to a coordinate axis in all ACR systems with two reactions and two species.

Each static and dynamic ACR motif (with two reactions and two species) has the following properties: the reactant polytope is parallel to a coordinate axis, the two reaction arrows are mutually parallel and they point towards each other. The three motifs are then characterized by the slope of the reaction vectors, which can be positive, zero, or negative. We will refer to the motifs accordingly as “positive/zero/negative slope motif”.

For each motif, we choose a particular network and show that when considered as an open system by allowing inflows, the ACR property is preserved in many cases. One representative of negative slope motif, see (3.2), is often studied as a simple example of a conservative system with static ACR. We show that not only does it have dynamic ACR, it has many other striking robustness properties.

4.2 Positive slope motif

Our first example is a network based on the motif . For this motif, we choose a particular reaction network whose true/non-flow reactions are $\{A+B\to 0,B\to A+2B\}$ , and we include the inflows $\{0\to A,0\to B\}$ with arbitrary time-dependent rates. So the reaction network is

	$\displaystyle A+B\xrightarrow{k_{1}}0,\quad$	$\displaystyle B\xrightarrow{k_{2}}A+2B,$
	$\displaystyle 0\xrightarrow{g_{a}(t)}A,\quad$	$\displaystyle 0\xrightarrow{g_{b}(t)}B.$

whose mass action system is

	$\displaystyle\dot{a}$	$\displaystyle=k_{1}b(k^{*}-a)+g_{a}(t),$		(4.1)
	$\displaystyle\dot{b}$	$\displaystyle=k_{1}b(k^{*}-a)+g_{b}(t),$		(4.1)

where $k^{*}=k_{2}/k_{1}$ .

Theorem 4.1.

Consider the mass action system (4.1) with time-dependent inflow rates $g_{a}:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ and $g_{b}:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ , such that $g_{a}(t)$ is bounded and

\int_{0}^{\infty}\left(g_{b}(t)-g_{a}(t)\right)dt=\infty,

(4.2)

then $a(t)\xrightarrow{t\to\infty}k^{*}$ .

Proof.

Since $\dot{b}-\dot{a}=g_{b}-g_{a}$ , by (4.2), we have that $b(t)-a(t)\xrightarrow{t\to\infty}\infty$ , which implies that $b(t)\xrightarrow{t\to\infty}\infty$ . Then the conditions of Theorem 3.8 hold for $a$ and from

\lim_{t\to\infty}\frac{g_{a}(t)}{b(t)}=0,

we conclude that $a(t)\xrightarrow{t\to\infty}k^{*}$ . ∎

Corollary 4.2.

Consider the mass action system (4.1) with constant inflows $g_{b}>g_{a}$ . Then $a(t)\xrightarrow{t\to\infty}k^{*}$ .

4.3 Negative slope motif

A representative of the motif was used to illustrate static ACR by Shinar and Feinberg [3]. We show here that this network has even stronger robustness properties than previously known. In particular, even in face of inflows which send the total concentration to infinity, the species with ACR still converges to a finite value, and moreover, the finite value is the same as its ACR value when there are no inflows.

Consider the following open reaction network:

	$\displaystyle A+B\xrightarrow{k_{1}}2B,$	$\displaystyle B\xrightarrow{k_{2}}A,$		(4.3)
	$\displaystyle 0\xrightarrow{g_{a}(t)}A,$	$\displaystyle 0\xrightarrow{g_{b}(t)}B,$		(4.3)

whose mass action system is

	$\displaystyle\dot{a}$	$\displaystyle=k_{1}b(k^{*}-a)+g_{a}(t),$		(4.4)
	$\displaystyle\dot{b}$	$\displaystyle=-k_{1}b(k^{*}-a)+g_{b}(t),$		(4.4)

where $k^{*}=k_{2}/k_{1}$ .

If $g_{a}$ and $g_{b}$ are constant, then the ACR value survives. We give a more general result below, allowing $g_{a}$ and $g_{b}$ to be arbitrary functions of time. One implication of the result is that the ACR value survives for arbitrary functions, provided that $g_{a}(t)$ does not increase too fast.

Theorem 4.3.

Consider the mass action system (4.4) with time-dependent inflow rates $g_{a}:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ and $g_{b}:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ , such that

G(t)\coloneqq\int_{0}^{t}\left(g_{a}(s)+g_{b}(s)\right)ds\xrightarrow{t\to\infty}\infty,

(4.5)

and

g_{a}(t)/G(t)\xrightarrow{t\to\infty}\alpha.

(4.6)

then $a(t)\xrightarrow{t\to\infty}k^{*}+\alpha/k_{1}$ for any $(a(0),b(0))\in\mathbb{R}^{2}_{>0}$ .

Proof.

We will first show that $b/G\to 1$ . Note that $\dot{a}+\dot{b}=g_{a}+g_{b}$ implies that $a(t)+b(t)=a(0)+b(0)+G(t),$ and so $a+b\xrightarrow{t\to\infty}\infty$ .

Define the following invertible change of coordinates:

	$\displaystyle\mathbb{R}_{\geq 0}^{2}\setminus\{(0,0)\}$	$\displaystyle\to\mathbb{R}_{>0}\times[0,1]$		(4.7)
	$\displaystyle(a,b)$	$\displaystyle\mapsto\left(x,\beta\right)=\left(a+b,b/(a+b)\right)$		(4.7)

In $(x,\beta)$ coordinates, the dynamical system (4.4) is:

	$\displaystyle\dot{x}$	$\displaystyle=g_{a}(t)+g_{b}(t),$		(4.8)
	$\displaystyle\dot{\beta}$	$\displaystyle=k_{1}\beta\left(x(1-\beta)-k^{*}\right)+\frac{(1-\beta)g_{b}-\beta g_{a}}{x}.$		(4.8)

The first equation has the solution $x(t)=x(0)+G(t)$ , which implies that $x(t)$ grows monotonically to infinity. It then follows that for any $\delta\in(0,\frac{1}{2})$ there is a time $T_{\delta}>0$ such that for all $t>T_{\delta}$ we have $\dot{\beta}(t)>1$ whenever $\beta(t)\in(\delta,1-\delta)$ , because under such assumptions $x(t)$ is large enough to imply that the positive term

k_{1}\beta\left(x(1-\beta)\right)

is much larger than all the negative terms on the right-hand side of the equation for $\dot{\beta}(t)$ . Moreover, we obtain that if $\beta(0)\in(\delta,1-\delta)$ then the solution $\beta(t)$ becomes larger than $1-\delta$ for all $t>T_{\delta}+1$ , which implies that $\beta(t)\to 1$ , and therefore $b/G\to 1$ .

Finally, we check the hypotheses of Theorem 3.8. We may rewrite the equation for $\dot{a}$ in (4.4) as

\displaystyle\dot{a}

\displaystyle=k_{1}b\left(\frac{\alpha}{k_{1}}+k^{*}-a\right)+G(t)\left(\frac{g_{a}(t)}{G(t)}-\alpha\frac{b}{G(t)}\right).

(4.9)

Clearly $b\to\infty$ and so $\int b\to\infty$ . Moreover, we have that

\lim_{t\to\infty}\frac{G(t)}{b(t)}\left(\frac{g_{a}(t)}{G(t)}-\alpha\frac{b}{G(t)}\right)=\lim_{t\to\infty}\left(\frac{g_{a}(t)}{G(t)}-\alpha\right)=0.

Therefore we may conclude from Theorem 3.8 that $a(t)\xrightarrow{t\to\infty}k^{*}+\alpha/k_{1}$ for any $(a(0),b(0))\in\mathbb{R}^{2}_{>0}$ . ∎

Remark 4.4.

Theorem 4.3 is formally about a two-dimensional non-autonomous system. We show in examples below that the theorem can be applied to study the dynamics of higher dimensional systems. See Examples 4.6, 4.7 and 4.19.

Corollary 4.5.

Consider the mass action system (4.4) with time-dependent inflow rates $g_{a}:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ and $g_{b}:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ , such that $g_{a}(t)$ is a polynomial of $t$ (that takes only positive values). This includes the case of constant inflows. Then $a(t)\xrightarrow{t\to\infty}k^{*}$ for any $(a(0),b(0))\in\mathbb{R}^{2}_{>0}$ .

Proof.

We only need to check the hypotheses of Theorem 4.3. Clearly (4.5) holds. Moreover for (4.6), we have

\displaystyle 0\leq\frac{g_{a}(t)}{G(t)}\leq\frac{g_{a}(t)}{\int_{0}^{t}g_{a}(s)ds}\xrightarrow{t\to\infty}0,

and so $\alpha=0$ . It follows that $a\xrightarrow{t\to\infty}k^{*}$ . ∎

Example 4.6.

Consider the following reaction network for some positive integer $d$

$\displaystyle 0$	$\displaystyle\to X_{1},$	(4.10)
$\displaystyle X_{j}$	$\displaystyle\to X_{j}+X_{j+1},\quad(1\leq j\leq d-1)$
$\displaystyle X_{d}+Y$	$\displaystyle\xrightarrow{k_{1}}2Y,$
$\displaystyle Y$	$\displaystyle\xrightarrow{k_{2}}X_{d}.$

Note that $d=1$ is allowed and corresponds to constant inflow of $X_{d}=X_{1}$ . For $j\in\{1,\ldots,d-1\}$ , $X_{j}\sim t^{j}$ and so Corollary 4.5 implies that $X_{d}\to k_{2}/k_{1}$ . ∎

Now we consider the case where the inflows are exponentially growing in time.

Example 4.7.

Consider the following reaction network

$\displaystyle Z\xrightarrow{\alpha}2Z,$	$\displaystyle\quad Z\to Z+A$	(4.11)
$\displaystyle A+B$	$\displaystyle\xrightarrow{k_{1}}2B,$
$\displaystyle B$	$\displaystyle\xrightarrow{k_{2}}A.$

The $Z$ could represent an exponentially growing virus which produces a toxin $A$ . The differential equation for $Z$ is $\dot{z}=\alpha z$ which has the solution $z(t)=e^{\alpha t}$ (assuming that $z(0)=1$ ) and so the concentrations of $A$ and $B$ evolve according to:

	$\displaystyle\dot{a}$	$\displaystyle=k_{1}b(k^{*}-a)+e^{\alpha t},$		(4.12)
	$\displaystyle\dot{b}$	$\displaystyle=-k_{1}b(k^{*}-a).$		(4.12)

Now

G(t)=\int_{0}^{t}e^{\alpha s}ds=\frac{e^{\alpha t}-1}{\alpha},

which clearly goes to infinity and moreover $g_{a}(t)/G(t)\to\alpha$ . By Theorem 4.3, it then follows that

a(t)\xrightarrow{t\to\infty}k^{*}+\alpha/k_{1}.

A similar result holds even when both $A$ and $B$ have exponentially growing inflows. For instance, for the system

	$\displaystyle\dot{a}$	$\displaystyle=k_{1}b(k^{*}-a)+e^{\alpha t},$		(4.13)
	$\displaystyle\dot{b}$	$\displaystyle=-k_{1}b(k^{*}-a)+e^{\beta t},$		(4.13)

we have that

\displaystyle\lim_{t\to\infty}a(t)=\begin{cases}k^{*}+\alpha/k_{1}&\mbox{ if }\alpha>\beta\\ k^{*}+\alpha/(2k_{1})&\mbox{ if }\alpha=\beta\\ k^{*}&\mbox{ if }\alpha<\beta.\end{cases}

∎

Remark 4.8.

It is possible to find an explicit solution for the case of arbitrary time-dependent inflow in $A$ , $g_{a}(t)$ and no inflow in $B$ . We give such a solution here. One can see from the form of the solution, even though it is explicit, that the limiting value of $a(t)$ given in Theorem 4.3 is not obtained easily. Moreover, we emphasize the additional shortcoming that the explicit solution is for the restricted situation of no inflow in $B$ .

Consider the following system:

\displaystyle\dot{a}=-k_{1}ab+k_{2}b+g_{a}(t),\quad\dot{b}=k_{1}ab-k_{2}b.

(4.14)

Let $G_{a}(t)=\int_{0}^{t}g_{a}(s)ds$ , so $\dot{a}+\dot{b}=g(t)$ has the solution $a(t)+b(t)=a(0)+b(0)+G_{a}(t)$ . We can use this to write an ODE in $a$ only:

\displaystyle\dot{a}=-k_{1}(a-k^{*})(a(0)+b(0)+G_{a}(t)-a)+g_{a}(t),

(4.15)

where $k^{*}=k_{2}/k_{1}$ .

Theorem 4.9.

Let $g_{a}(t)$ be an arbitrary function of time $t$ and $G_{a}(t)=\int_{0}^{t}g_{a}(s)~{}ds$ . For any initial value $a(0)=a_{0}\geq 0$ and $b(0)=b_{0}\geq 0$ ,

\displaystyle\dot{a}=-k_{1}(a-k)(a_{0}+b_{0}+G_{a}(t)-a)+g_{a}(t)

(4.16)

has the solution

\displaystyle a(t)=a_{0}+b_{0}+G_{a}(t)-\frac{b_{0}q(t)}{1+k_{1}b_{0}Q(t)}

(4.17)

for all time $t\geq 0$ where

\displaystyle q(t)=\exp\left[k_{1}(a_{0}+b_{0}-k)t+k_{1}\int_{0}^{t}G_{a}(s)~{}ds\right],

(4.18)

and $Q(t)=\int_{0}^{t}q(s)~{}ds$ .

Proof.

The proof is a simple verification. From their definitions, $q(0)=1,G_{a}(0)=0,Q(0)=0$ and so

a_{0}+b_{0}+G_{a}(0)-\frac{b_{0}q(0)}{1+k_{1}b_{0}Q(0)}=a_{0}+b_{0}-b_{0}=a_{0},

so the initial condition is satisfied. We note that $\dot{q}(t)=k_{1}q(t)\left(a_{0}+b_{0}-k+G_{a}(t)\right)$ and so

\displaystyle\dot{a}(t)

\displaystyle=g_{a}(t)-b_{0}\frac{k_{1}q(t)\left(a_{0}+b_{0}-k+G_{a}(t)\right)}{1+k_{1}b_{0}Q(t)}+k_{1}\left(\frac{b_{0}q(t)}{1+k_{1}b_{0}Q(t)}\right)^{2},

and therefore

	$\displaystyle\dot{a}+k_{1}(a-k)(a_{0}+b_{0}+G_{a}(t)-a)-g_{a}(t)$
	$\displaystyle=\dot{a}+k_{1}\left(a_{0}+b_{0}+G_{a}(t)-\frac{b_{0}q(t)}{1+k_{1}b_{0}Q(t)}-k\right)\left(\frac{b_{0}q(t)}{1+k_{1}b_{0}Q(t)}\right)-g_{a}(t)=0,$

which completes the verification. ∎

4.4 Negative slope motif with inflows and outflows

We can extend the results in the previous section even to the situation where one or both species are in outflow. We consider the case where both $A$ and $B$ are in outflow and at equal rates. Outflows with equal rates for all species is a reasonable (and standard) assumption for chemostats. We will allow the inflows to be arbitrary functions of time. Suppose we have the following reaction network:

	$\displaystyle 0\xrightarrow{g_{a}(t)}A,\quad$	$\displaystyle 0\xrightarrow{g_{b}(t)}B,$
	$\displaystyle A\xrightarrow{\ell}0,\quad$	$\displaystyle B\xrightarrow{\ell}0,$
	$\displaystyle A+B\xrightarrow{k_{1}}2B,\quad$	$\displaystyle B\xrightarrow{k_{2}}A,$

whose mass action system is:

	$\displaystyle\dot{a}$	$\displaystyle=k_{1}b(k^{*}-a)+g_{a}(t)-\ell a,$		(4.19)
	$\displaystyle\dot{b}$	$\displaystyle=-k_{1}b(k^{*}-a)+g_{b}(t)-\ell b.$		(4.19)

Refer to caption — Figure 2: The reaction network $\{A+B\xrightarrow{k_{1}}2B,~{}B\xrightarrow{k_{2}}A\}$ coupled with an inflow $0\xrightarrow{g_{a}(t)}A$ (with $g_{a}(t)\xrightarrow{t\to\infty}\infty$ ) and outflows $A\xrightarrow{\ell}0,~{}B\xrightarrow{\ell}0$ . When the inflow $g_{a}(t)$ is sub-exponential, $a$ converges to a value that is equivalent to changing the rate constant of the reaction $B\to A$ from $k_{2}$ to $k_{2}+\ell$ . When the inflow $g_{a}(t)=\exp(\alpha t)$ is exponential, $a$ converges to a value that is equivalent to changing the rate constant of the reaction $B\to A$ from $k_{2}$ to $k_{2}+\alpha+\ell$ .

Theorem 4.10.

Consider (4.19) with $g_{a},g_{b}:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ such that $g(t)\coloneqq g_{a}(t)+g_{b}(t)\xrightarrow{t\to\infty}\infty$ , and $\ell>0$ . Assume that if $g_{b}\equiv 0$ then $b(0)>0$ . Let

H(t)\coloneqq e^{-\ell t}\int_{0}^{t}g(s)e^{\ell s}~{}ds.

If $\alpha\coloneqq\lim_{t\to\infty}g_{a}(t)/H(t)$ exists and is finite, then

a(t)\xrightarrow{t\to\infty}k^{*}+\frac{\alpha}{k_{1}}.

Proof.

We will first show that $b/H\to 1$ . Define the following invertible change of coordinates:

	$\displaystyle\mathbb{R}_{\geq 0}^{2}\setminus\{(0,0)\}$	$\displaystyle\to\mathbb{R}_{>0}\times[0,1]$
	$\displaystyle(a,b)$	$\displaystyle\mapsto\left(x,\beta\right)=\left(a+b,b/(a+b)\right).$

In $(x,\beta)$ coordinates, the dynamical system (4.19) is

	$\displaystyle\dot{x}$	$\displaystyle=g(t)-\ell x$
	$\displaystyle\dot{\beta}$	$\displaystyle=k_{1}\beta\left(x(1-\beta)-k^{*}\right)+\frac{(1-\beta)g_{b}-\beta g_{a}}{x}$

The equation for $\dot{x}$ can be solved explicitly using an integrating factor:

x(t)=e^{-\ell t}\left(\int_{0}^{t}g(s)e^{\ell s}~{}ds+x(0)\right)=x(0)e^{-\ell t}+H(t).

We now show that $x\to\infty$ . By hypothesis, we have that $g\xrightarrow{t\to\infty}\infty$ . So for any $N>0$ , there is a $T_{N}>(\ln 2)/\ell$ such that $g(t)>N$ for $t>T_{N}$ . So,

	$\displaystyle x(2T_{N})$	$\displaystyle\geq H(2T_{N})\geq e^{-2\ell T_{N}}\int_{T_{N}}^{2T_{N}}g(s)e^{\ell s}~{}ds$
		$\displaystyle\geq e^{-2\ell T_{N}}N\int_{T_{N}}^{2T_{N}}e^{\ell s}~{}ds=N\frac{1-\exp(-\ell T_{N})}{\ell}\geq\frac{N}{2\ell}.$

Since $N$ is arbitrary, $x\xrightarrow{t\to\infty}\infty$ .

Now note that the equation for $\dot{\beta}$ is the same as the corresponding equation in Theorem 4.3 (in particular, the outflow parameter $\ell$ does not appear in $\dot{\beta}$ ). Therefore, the same reasoning as in Theorem 4.3 gives us that $b/x\xrightarrow{t\to\infty}1$ and $b/H\xrightarrow{t\to\infty}1$ . Moreover, $a/x\to 0$ and

\frac{a}{b}=\frac{a/x}{b/x}\to 0.

Finally, we check the hypotheses of Theorem 3.8. Clearly $b\to\infty$ and so $\int b\to\infty$ . Moreover, we have that

\lim_{t\to\infty}\frac{-\ell a+g_{a}(t)}{k_{1}b(t)}=\lim_{t\to\infty}\frac{1}{k_{1}}\left(-\ell\frac{a}{b}+\frac{g_{a}(t)/H(t)}{b/H(t)}\right)=\frac{\alpha}{k_{1}}.

Therefore, we conclude from Theorem 3.8 that $a(t)\xrightarrow{t\to\infty}k^{*}+\alpha/k_{1}$ for any $b(0)>0$ or $g_{b}\not\equiv 0$ . ∎

Remark 4.11.

In Theorem 4.10, we assumed equal outflows of species $A$ and $B$ , $\ell_{a}=\ell_{b}=\ell$ , a standard assumption for a chemostat. Simulations show that even when this assumption is relaxed and the outflow rates are allowed to be different, the property of dynamic ACR in $A$ is preserved. We conjecture that the ACR value is obtained by replacing $\ell$ with $\ell_{b}$ , the outflow of species $B$ , in the expression for $H(t)$ in Theorem 4.10.

Corollary 4.12.

Assume the hypotheses and notation of Theorem 4.10. Suppose that $g_{a}$ is differentiable. Then

\alpha=\lim_{t\to\infty}\left[\ell\frac{g_{a}(t)}{g(t)}+\frac{g_{a}^{\prime}(t)}{g(t)}\right],

and so

a(t)\xrightarrow{t\to\infty}k^{*}+\frac{1}{k_{1}}\lim_{t\to\infty}\left[\ell\frac{g_{a}(t)}{g(t)}+\frac{g_{a}^{\prime}(t)}{g(t)}\right].

Proof.

Note that

\displaystyle\alpha=\lim_{t\to\infty}\frac{g_{a}(t)}{H(t)}=\lim_{t\to\infty}\frac{g_{a}(t)e^{\ell t}}{\int_{0}^{t}g(s)e^{\ell s}ds}.

For any $g_{a}(t)\not\equiv 0$ , both the numerator and the denominator go to infinity, so the result then follows from applying L’Hospital’s rule. ∎

Example 4.13.

Consider the reaction network (4.10) in Example 4.6 with $d\geq 2$ (so that $g_{a}(t)\to\infty$ ), and additionally include outflows $X_{d}\xrightarrow{\ell}0,~{}Y\xrightarrow{\ell}0,$ so the complete reaction network is

$\displaystyle 0$	$\displaystyle\to X_{1},$	(4.20)
$\displaystyle X_{j}$	$\displaystyle\to X_{j}+X_{j+1},\quad(1\leq j\leq d-1)$
$\displaystyle X_{d}+Y$	$\displaystyle\xrightarrow{k_{1}}2Y,$
$\displaystyle Y$	$\displaystyle\xrightarrow{k_{2}}X_{d},$
$\displaystyle X_{d}$	$\displaystyle\xrightarrow{\ell}0,\quad Y\xrightarrow{\ell}0.$

As earlier, for $j\in\{1,\ldots,d-1\}$ , $X_{j}\sim t^{j}$ . Choose rate constants and initial values so that $X_{j}=t^{j}$ for $j\in\{1,\ldots,d-1\}$ . So the concentrations for $X_{d}$ and $Y$ follow:

	$\displaystyle\dot{x}_{d}$	$\displaystyle=k_{1}y(k^{*}-x_{d})+t^{d-1}-\ell x_{d},$
	$\displaystyle\dot{y}$	$\displaystyle=-k_{1}y(k^{*}-x_{d})-\ell y.$

In the notation of Theorem 4.10, $g_{x_{d}}(t)=g(t)=t^{d-1}$ and by Corollary 4.12

\displaystyle\alpha=\ell+\lim_{t\to\infty}\frac{g_{a}^{\prime}(t)}{g(t)}=\ell,

and so

x_{d}\xrightarrow{t\to\infty}k^{*}+\frac{\ell}{k_{1}}.

See left panel of Figure 2 for a representative trajectory. ∎

Example 4.14.

Consider the reaction network (4.11) in Example 4.7 with outflows included, so the complete reaction network is

$\displaystyle Z\xrightarrow{\alpha}2Z,$	$\displaystyle\quad Z\to Z+A$	(4.21)
$\displaystyle A+B$	$\displaystyle\xrightarrow{k_{1}}2B,$
$\displaystyle B$	$\displaystyle\xrightarrow{k_{2}}A,$
$\displaystyle B\xrightarrow{\ell}0,$	$\displaystyle\quad A\xrightarrow{\ell}0.$

The concentrations of $A$ and $B$ evolve according to:

	$\displaystyle\dot{a}$	$\displaystyle=k_{1}b(k^{*}-a)+e^{\alpha t}-\ell a,$		(4.22)
	$\displaystyle\dot{b}$	$\displaystyle=-k_{1}b(k^{*}-a)-\ell b.$		(4.22)

Noting that $g_{a}(t)=g(t)=e^{\alpha t}$ , by Corollary 4.12,

a\xrightarrow{t\to\infty}k^{*}+\frac{\ell}{k_{1}}+\frac{\alpha}{k_{1}}.

See right panel of Figure 2 for a representative trajectory. ∎

4.5 A representative of the negative slope motif with stronger robustness properties

Another representative of the negative slope motif has even stronger robustness properties. Consider the following open reaction network:

	$\displaystyle A+2B\xrightarrow{k_{1}}3B,$	$\displaystyle 2B\xrightarrow{k_{2}}A+B,$		(4.23)
	$\displaystyle 0\xrightarrow{g_{a}(t)}A,$	$\displaystyle 0\xrightarrow{g_{b}(t)}B,$		(4.23)

whose mass action system is

	$\displaystyle\dot{a}$	$\displaystyle=k_{1}b^{2}(k^{*}-a)+g_{a}(t),$		(4.24)
	$\displaystyle\dot{b}$	$\displaystyle=-k_{1}b^{2}(k^{*}-a)+g_{b}(t),$		(4.24)

where $k^{*}=k_{2}/k_{1}$ .

Theorem 4.15.

Consider the mass action system (4.24) with time-dependent inflow rates $g_{a}:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ and $g_{b}:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ , such that

G(t)\coloneqq\int_{0}^{t}\left(g_{a}(s)+g_{b}(s)\right)ds\xrightarrow{t\to\infty}\infty,

(4.25)

and

g_{a}(t)/(G(t))^{2}\xrightarrow{t\to\infty}\alpha.

(4.26)

Then $a(t)\xrightarrow{t\to\infty}k^{*}+\alpha/k_{1}$ for any $(a(0),b(0))\in\mathbb{R}^{2}_{>0}$ .

Proof.

The proof is analogous to that of Theorem 4.3. We will first show that $b/G\to 1$ . Note that $\dot{a}+\dot{b}=g_{a}+g_{b}$ implies that $a(t)+b(t)=a(0)+b(0)+G(t),$ and so $a+b\xrightarrow{t\to\infty}\infty$ . After changing to $(x,\beta)$ coordinates defined in (4.7), the dynamical system (4.4) is:

	$\displaystyle\dot{x}$	$\displaystyle=g_{a}(t)+g_{b}(t),$		(4.27)
	$\displaystyle\dot{\beta}$	$\displaystyle=k_{1}\beta^{2}x\left(x(1-\beta)-k^{*}\right)+\frac{(1-\beta)g_{b}-\beta g_{a}}{x}.$		(4.27)

k_{1}\beta^{2}x\left(x(1-\beta)\right)

Finally, we check the hypotheses of Theorem 3.8. We may rewrite the equation for $\dot{a}$ in (4.4) as

\displaystyle\dot{a}

\displaystyle=k_{1}b^{2}\left(\frac{\alpha}{k_{1}}+k^{*}-a\right)+G^{2}(t)\left(\frac{g_{a}(t)}{G^{2}(t)}-\alpha\frac{b^{2}}{G^{2}(t)}\right).

(4.28)

Clearly $b\to\infty$ and so $\int b\to\infty$ . Moreover, we have that

\lim_{t\to\infty}\frac{G^{2}(t)}{b^{2}(t)}\left(\frac{g_{a}(t)}{G^{2}(t)}-\alpha\frac{b^{2}}{G^{2}(t)}\right)=\lim_{t\to\infty}\left(\frac{g_{a}(t)}{G^{2}(t)}-\alpha\right)=0.

Therefore we may conclude from Theorem 3.8 that $a(t)\xrightarrow{t\to\infty}k^{*}+\alpha/k_{1}$ for any $(a(0),b(0))\in\mathbb{R}^{2}_{>0}$ . ∎

Theorem 4.16.

Consider the mass action system (4.24) with time-dependent inflow rates $g_{a}:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ and $g_{b}:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ , such that

G(t)\coloneqq\int_{0}^{t}\left(g_{a}(s)+g_{b}(s)\right)ds\xrightarrow{t\to\infty}\infty.

(4.29)

We have that $a(t)\xrightarrow{t\to\infty}k^{*}$ for any $(a(0),b(0))\in\mathbb{R}^{2}_{>0}$ if one of the following holds:

1.

$g_{a}$ is bounded,
2.

$g_{a}(t)\to\infty$ and $\displaystyle\frac{\dot{g}_{a}(t)}{g_{a}(t)^{1/2}g(t)}\to 0$ .

Proof.

If $g_{a}$ is bounded, then $g_{a}(t)/(G(t))^{2}\xrightarrow{t\to\infty}0$ and so the result follows immediately from Theorem 4.15. Suppose that $g_{a}(t)\xrightarrow{t\to\infty}\infty$ and $\displaystyle\frac{\dot{g}_{a}(t)}{g_{a}(t)^{1/2}g(t)}\to 0$ . Since $g_{a}(t)\to\infty$ , clearly $G(t)\to\infty$ and so by applying L’Hospital’s rule, we have that

\displaystyle\lim_{t\to\infty}\frac{\sqrt{g_{a}(t)}}{G(t)}=\lim_{t\to\infty}\frac{\frac{d}{dt}\sqrt{g_{a}(t)}}{\frac{d}{dt}G(t)}=\lim_{t\to\infty}\frac{\dot{g}_{a}(t)}{2g_{a}(t)^{1/2}g(t)}=0,

and so $g_{a}(t)/(G(t))^{2}\to 0$ . Then from Theorem 4.15, $a(t)\xrightarrow{t\to\infty}k^{*}$ . ∎

Remark 4.17.

When $g_{a}(t)\xrightarrow{t\to\infty}\infty$ , the limit of $\displaystyle\frac{\dot{g}_{a}(t)}{g_{a}(t)^{1/2}g(t)}$ cannot be positive. Indeed suppose that $\displaystyle\frac{\dot{g}_{a}(t)}{g_{a}(t)^{1/2}g(t)}\xrightarrow{t\to\infty}\alpha>0$ , then after some finite time $T_{c}>0$ , we have that for some positive constant $c>0$ ,

\frac{d}{dt}g_{a}(t)>c\alpha g_{a}(t)^{1/4}g(t)\geq c\alpha g_{a}(t)^{5/4},

and the differential equation $\dot{g}_{a}(t)=c\alpha g_{a}(t)^{5/4}$ has a finite time blow up.

Remark 4.18.

A simpler result is obtained by replacing the term $\displaystyle\frac{\dot{g}_{a}(t)}{g_{a}(t)^{1/2}g(t)}$ in Theorem 4.16 and Remark 4.17 with the term $\displaystyle\frac{\dot{g}_{a}(t)}{g_{a}(t)^{3/2}}$ .

Example 4.19.

Consider the following reaction network for some positive integer $m$

$\displaystyle Z_{1}$	$\displaystyle\xrightarrow{1}2Z_{1},$	(4.30)
$\displaystyle Z_{2}+Z_{1}$	$\displaystyle\xrightarrow{1}2Z_{2}+Z_{1},$
$\displaystyle Z_{2}$	$\displaystyle\xrightarrow{1}Z_{2}+A$
$\displaystyle A+2B$	$\displaystyle\xrightarrow{k_{1}}3B,$
$\displaystyle 2B$	$\displaystyle\xrightarrow{k_{2}}A+B.$

We assumed that the rate constants for some of the equations are $1$ for simplicity. The variable $z_{1}$ satisfies the differential equation $\dot{z}_{1}=z_{1}$ , which has an exponentially growing solution, $z_{1}=e^{t}$ , assuming initially $z_{1}(0)=1$ . Then $z_{2}$ satisfies the differential equation

\dot{z}_{2}=z_{1}z_{2},

which has the solution:

z_{2}(t)=e^{e^{t}},

assuming that $z_{2}(0)=e$ . The ODE system satisfied by the species $A$ and $B$ is:

	$\displaystyle\dot{a}$	$\displaystyle=k_{1}b^{2}(k^{*}-a)+e^{e^{t}},$		(4.31)
	$\displaystyle\dot{b}$	$\displaystyle=-k_{1}b^{2}(k^{*}-a).$		(4.31)

Let $g_{a}(t)\coloneqq e^{e^{t}}$ , then clearly $g_{a}(t)\to\infty$ and

\frac{\dot{g}_{a}(t)}{g_{a}(t)^{3/2}}=\frac{e^{t}e^{e^{t}}}{e^{3e^{t}/2}}=e^{t-e^{t}/2}\xrightarrow{t\to\infty}0.

Then, by Theorem 4.16, $a\to k^{*}$ . ∎

4.6 Zero slope motif

We consider the motif and the following open reaction network:

	$\displaystyle B\mathrel{\mathop{\vbox{\offinterlineskip\halign{\hfil#\hfil\cr\hphantom{$\scriptstyle\mspace{8.0mu}{k_{2}}\mspace{8.0mu}$}\cr\rightarrowfill\cr\vrule height=0.0pt,width=20.00003pt\cr\leftarrowfill\cr\hphantom{$\scriptstyle\mspace{8.0mu}{k_{1}}\mspace{8.0mu}$}\cr\kern-1.29167pt\cr}}}\limits^{k_{2}}_{k_{1}}}A+B,$			(4.32)
	$\displaystyle 0\xrightarrow{g_{a}(t)}A,$	$\displaystyle 0\xrightarrow{g_{b}(t)}B.$		(4.32)

whose mass action system is

	$\displaystyle\dot{a}$	$\displaystyle=k_{1}b(k^{*}-a)+g_{a}(t),$		(4.33)
	$\displaystyle\dot{b}$	$\displaystyle=g_{b}(t),$		(4.33)

where $k^{*}=k_{2}/k_{1}$ .

Theorem 4.20.

Consider (4.33) with $g_{a},g_{b}:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ such that $\int_{0}^{t}b(s)ds\xrightarrow{t\to\infty}\infty$ . Then

a(t)\xrightarrow{t\to\infty}k^{*}+\frac{1}{k_{1}}\lim_{t\to\infty}\frac{g_{a}(t)}{b(0)+\int_{0}^{t}g_{b}(s)ds},

provided the limit exists as a finite value.

Proof.

The result follows immediately from applying Theorem 3.8. ∎

Remark 4.21.

It is worth contrasting the dynamic ACR properties of the negative slope motif with the other motifs. The negative slope motif has particularly striking robustness properties when inflows are added, which is not the case for the other motifs. This is further argument in favor of the attention that the negative slope motif has attracted, besides having a positive conservation law and being relatively simple in nature.

4.7 Enzyme catalysis with bifunctional enzyme

Bifunctional enzymes are those which have two different (usually opposing) functions, for instance helping both to catalyze production as well as elimination of a substrate [8, 21]. For instance, a bifunctional enzyme might facilitate phosphorylation and when in a different form, facilitate dephosphorylation. Bifunctional enzymes have been shown to play a critical role in generating static absolute concentration robustness [9, 21]. Here we show that an enzyme catalysis network with a bifunctional enzyme can produce dynamic absolute concentration robustness as well.

We consider a mass action system with inflows that results from the reaction network depicted below.

$\displaystyle 0\xrightarrow{g_{x}}X,~{}0\xrightarrow{g_{y}}Y,{}$	$\displaystyle 0\xrightarrow{g_{e}}E,~{}0\xrightarrow{g_{c}}C,$	(4.34)
$\displaystyle X+E\stackrel{{\scriptstyle[}}{{k}}_{2}]{k_{1}}{\rightleftarrows}$	$\displaystyle C\stackrel{{\scriptstyle k_{3}}}{{\to}}Y+E,$
$\displaystyle Y+C$	$\displaystyle\xrightarrow{k_{4}}X+C.$

The system of ODEs for the mass action system is

$\displaystyle\dot{x}$	$\displaystyle=-k_{1}xe+k_{2}c+k_{4}cy+g_{x},$	(4.35)
$\displaystyle\dot{y}$	$\displaystyle=k_{3}c-k_{4}cy+g_{y},$
$\displaystyle\dot{e}$	$\displaystyle=-k_{1}xe+(k_{2}+k_{3})c+g_{e},$
$\displaystyle\dot{c}$	$\displaystyle=k_{1}xe-(k_{2}+k_{3})c+g_{c}.$

We note that the variable $y$ is in power-engine-load form:

\displaystyle\dot{y}={\color[rgb]{0.5,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.5,0.5,0}\pgfsys@color@cmyk@stroke{0}{0}{1}{.5}\pgfsys@color@cmyk@fill{0}{0}{1}{.5}k_{4}c}{\color[rgb]{1,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,1}\pgfsys@color@cmyk@stroke{0}{1}{0}{0}\pgfsys@color@cmyk@fill{0}{1}{0}{0}\left(\frac{k_{3}}{k_{4}}-y\right)}+{\color[rgb]{0,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,1}\pgfsys@color@cmyk@stroke{1}{0}{0}{0}\pgfsys@color@cmyk@fill{1}{0}{0}{0}g_{y}},

and therefore is a candidate for dynamic ACR. The key to proving dynamic ACR for different inflows or no inflows is to show that the “power” variable $c$ is eventually bounded away from $0$ for any initial value. We will do this here for a variety of cases with inflows. When there is no inflow in $Y$ ( $g_{y}=0$ ) and assuming that $Y$ is a dynamic ACR species, then it is easy to see that the dynamic ACR value must be $k^{*}\coloneqq k_{3}/k_{4}$ , the same as its static ACR value. In the case of no inflows (in any species $g_{x}=g_{y}=g_{e}=g_{c}=0$ ), one can further show that $y$ has dynamic ACR, that is for any positive initial value with $T_{x}\coloneqq x(0)+y(0)+c(0)\geq k^{*}$ , we have that $y\to k^{*}$ . In other words, every initial value compatible with the hyperplane $\{y=k^{*}\}$ results in convergence to the hyperplane. We defer the proof of this result to future work because it requires additional ideas not directly related to the main theme of the paper. We focus here on certain interesting inflow cases.

In the theorem below, we show that $y$ can converge to the same value $k^{*}$ even in cases when some (or all) of the inflows are present.

Theorem 4.22.

Consider (4.35) with $g_{c}>0$ (while the other inflow rates $g_{x},g_{y},g_{e}$ are arbitrary nonnegative constants). Then $y\xrightarrow{t\to\infty}k^{*}\coloneqq k_{3}/k_{4}$ .

Proof.

Introduce the invertible change of variables $(x,y,e,c)\mapsto(y,c,T_{x},T_{e})\in\mathbb{R}^{4}_{\geq 0}$ , where $T_{x}\coloneqq x+y+c$ and $T_{e}\coloneqq e+c$ . The new variables satisfy

$\displaystyle\dot{y}$	$\displaystyle=k_{4}c\left(k^{*}-y\right)+g_{y}$	(4.36)
$\displaystyle\dot{c}$	$\displaystyle=k_{1}\left((T_{x}-y-c)(T_{e}-c)-kc\right)+g_{c}$
$\displaystyle\dot{T}_{x}$	$\displaystyle=g_{x}+g_{y}+g_{c}=:g_{1}$
$\displaystyle\dot{T}_{e}$	$\displaystyle=g_{e}+g_{c}=:g_{2}$

where $k\coloneqq(k_{2}+k_{3})/k_{1}$ , $k^{*}\coloneqq k_{3}/k_{4}$ , as well as $0\leq c(t)\leq T_{e}(t)$ and $0\leq c(t)+y(t)\leq T_{x}(t)$ for all $t\geq 0$ . The hypothesis $g_{c}>0$ implies that $g_{1}>0$ and $g_{2}>0$ . The evolution of $T_{x}$ and $T_{e}$ is uncoupled from that of $y$ and $c$ , and is easily seen to have the solution $T_{x}(t)=T_{x}(0)+g_{1}t$ and $T_{e}(t)=T_{e}(0)+g_{2}t$ . So it is natural to think of (4.36) as a two-dimensional non-autonomous system in the variables $(y,c)$ .

Since the first term of $\dot{c}$ is nonnegative, we have that $\dot{c}\geq-(k_{2}+k_{3})c+g_{c}$ which is positive for any $c<g_{c}/(k_{2}+k_{3})$ and so after some finite time, $c(t)>g_{c}/(k_{2}+k_{3})$ . But this implies that $\dot{y}<0$ for any sufficiently large $y$ , and so $y$ is bounded above by some $\bar{y}>0$ after some finite time.

Now suppose that

c(t)<\frac{1}{2}\left(\min\left\{T_{x}(t)-\bar{y},T_{e}(t)\right\}-\frac{k}{2}\right).

Then

	$\displaystyle\dot{c}$	$\displaystyle=k_{1}\left((T_{x}(t)-y-c)(T_{e}(t)-c)-kc\right)+g_{c}$
		$\displaystyle>k_{1}\left((T_{x}(t)-\bar{y}-c)(T_{e}(t)-c)-kc\right)+g_{c}$
		$\displaystyle>k_{1}\left(\left(c+k/2\right)^{2}-kc\right)+g_{c}$
		$\displaystyle=k_{1}\left(c^{2}+k^{2}/4\right)+g_{c},$

which implies that $c\to\infty$ since $T_{x}(t)-\bar{y}\to\infty$ and $T_{e}(t)\to\infty$ .

The result $y\to k^{*}$ then follows from applying Theorem 3.8. ∎

We now present a case where specific inflows result in $y$ converging to a finite value, which is however not the same as $k^{*}$ (unless $g_{y}=0$ ), the ACR value of $y$ when there are no inflows.

Theorem 4.23.

Consider (4.35) with $g_{c}=0$ , $g_{e}=0$ , $g_{y}\geq 0$ and $g_{x}+g_{y}>0$ . Then

y\xrightarrow{t\to\infty}k^{*}+\frac{g_{y}}{k_{4}(e(0)+c(0))}.

In particular, $y\to k^{*}$ if $g_{y}=0$ .

Proof.

Introduce the invertible change of variables $(x,y,e,c)\mapsto(y,c,T_{x},T_{e})\in\mathbb{R}^{4}_{\geq 0}$ , where $T_{x}\coloneqq x+y+c$ and $T_{e}\coloneqq e+c$ . The new variables satisfy

$\displaystyle\dot{y}$	$\displaystyle=k_{4}c\left(k^{*}-y\right)+g_{y}$	(4.37)
$\displaystyle\dot{c}$	$\displaystyle=k_{1}\left((T_{x}(t)-y-c)(T_{e}-c)-kc\right)$
$\displaystyle\dot{T}_{x}$	$\displaystyle=g_{x}+g_{y}=:g_{1}>0$
$\displaystyle T_{e}(t)$	$\displaystyle=T_{e}(0)\coloneqq T_{e}$

where $k\coloneqq(k_{2}+k_{3})/k_{1}$ , $k^{*}\coloneqq k_{3}/k_{4}$ , as well as $0\leq c(t)\leq T_{e}$ and $0\leq c(t)+y(t)\leq T_{x}(t)$ for all $t\geq 0$ .

Suppose first that $y\leq k^{*}$ for all time. For arbitrary $\varepsilon$ , suppose that $c<T_{e}-\varepsilon$ . Then we have that

	$\displaystyle\dot{c}$	$\displaystyle=k_{1}\left((T_{x}(t)-y-c)(T_{e}-c)-kc\right)$
		$\displaystyle>k_{1}\left((T_{x}(t)-k^{*}-c)\varepsilon-kc\right)$
		$\displaystyle>k_{1}\left((T_{x}(t)-k^{*}-T_{e})\varepsilon-kT_{e}\right).$

Since $T_{x}(t)\uparrow\infty$ , we have that the right hand side is arbitrarily large after some finite time, which implies that $c(t)\to T_{e}$ .

Now suppose that $y(t)>k^{*}$ for some $t$ . It is clear that $\{y>k^{*}\}$ is invariant for any $g_{y}\geq 0$ . We may therefore assume, without loss of generality, that $y(t)>k^{*}$ for all $t\geq 0$ . For convenience we add $0=0\cdot c=(k_{3}-k_{4}k^{*})c$ to $\dot{x}$ and rewrite the original system (4.35) as

$\displaystyle\dot{x}$	$\displaystyle=-k_{1}xe+(k_{2}+k_{3})c-k_{4}c(k^{*}-y)+g_{x},$	(4.38)
$\displaystyle\dot{y}$	$\displaystyle=k_{4}c(k^{*}-y)+g_{y},$
$\displaystyle\dot{e}$	$\displaystyle=-k_{1}xe+(k_{2}+k_{3})c,$
$\displaystyle\dot{c}$	$\displaystyle=k_{1}xe-(k_{2}+k_{3})c.$

Let $u(t)\coloneqq k_{1}xe-(k_{2}+k_{3})c=k_{1}(xe-kc),~{}v(t)\coloneqq-k_{4}c(k^{*}-y)$ , and so the previous system in terms of $u$ and $v$ is:

\displaystyle\dot{x}=-u(t)+v(t)+g_{x},\quad\dot{y}=-v(t)+g_{y},\quad\dot{e}=-u(t),\quad\dot{c}=u(t).

(4.39)

Note that $y(t)>k^{*}$ is equivalent to $v(t)>0$ . Furthermore,

$\displaystyle\dot{u}(t)$	$\displaystyle=k_{1}(\dot{x}e+x\dot{e}-k\dot{c})$	(4.40)
	$\displaystyle=k_{1}(-eu+ev+g_{x}e-xu-ku)=k_{1}(-u(e+x+k)+e(v+g_{x}))$
	$\displaystyle>-k_{1}(e+x+k)u$

This shows that either $u(t)>0$ for all time after some finite time has passed, or $u(t)\uparrow 0$ . We argue that in both cases $c$ must converge to a positive limit.

First suppose that $u(t)\uparrow 0$ . Then $\dot{x}\geq 0$ and $\dot{e}\geq 0$ , and so $x(t)e(t)$ is nondecreasing for any $t\geq 0$ . If $xe\uparrow\infty$ , then $u(t)=k_{1}(xe-kc)\uparrow\infty$ since $c$ is bounded above. But this contradicts $u(t)=k_{1}(xe-kc)\uparrow 0$ . So $xe$ is bounded above and therefore converges to a positive limit that is greater than $x(0)e(0)$ . But then from the equation $\dot{c}=k_{1}xe-(k_{2}+k_{3})c$ , we have that $c(t)$ converges to a positive limit that is greater than $x(0)e(0)/k$ .

Now suppose that $u(t)>0$ for all large enough $t$ . Then $\dot{c}=u(t)>0$ for large enough $t$ , so that $c(t)$ is increasing. Since $c(t)$ is bounded above by $T_{e}$ , $c(t)$ must converge to a positive limit.

We conclude that $c$ converges to a positive limit in every case. From the equation $\dot{y}=k_{4}c(k^{*}-y)+g_{y}$ , we conclude that $y$ converges to a positive limit. If $c$ converges to a limit less than $T_{e}$ , then from $\dot{c}=k_{1}\left((T_{x}(t)-y-c)(T_{e}-c)-kc\right)$ , $c$ increases without bound because $T_{x}$ increases without bound. This is a contradiction. Therefore $c\to T_{e}$ . We may therefore conclude from Theorem 3.8 that

y\xrightarrow{t\to\infty}k^{*}+\frac{g_{y}}{k_{4}(e(0)+c(0))}.

∎

5 Discussion

The property of dynamic ACR depends on global dynamics and this makes it difficult to prove it even in simple systems. It is helpful to have either network conditions or ODE characterizations. Extending our previous work on network motifs which was only in two dimensions, we provide a power-engine-load form applicable in arbitrary dimensions. We give convergence results and apply these to certain network models in the chemostat setting, where inflows and sometimes outflows are present. Even when the overall trajectory does not converge, we show that the ACR variable converges to a definite finite quantity, the so-called ACR value. Moreover, in many cases this is the same value that the system would converge to when taken as a closed system, i.e. without any inflows and outflows.

This work has implications for understanding the dynamics of systems which are far from equilibrium. We focus here on showing dynamic ACR in a few representative examples. In future work, we apply some of these results to analyze more complex and biochemically realistic systems.

Acknowledgments

B.J. was supported by NSF grant DMS-2051498. G.C. was supported by NSF grant DMS-2051568 and by a Simons Foundation fellowship. We thank the referee for careful reading and helpful comments.

References

[1] Badal Joshi and Gheorghe Craciun. Foundations of static and dynamic absolute concentration robustness. Journal of Mathematical Biology, 85(53), 2022.
[2] Badal Joshi and Gheorghe Craciun. Reaction network motifs for static and dynamic absolute concentration robustness. SIAM Journal on Applied Dynamical Systems, 22(2):501–526, 2023.
[3] Guy Shinar and Martin Feinberg. Structural sources of robustness in biochemical reaction networks. Science, 327(5971):1389–1391, 2010.
[4] Frank D Russo and Thomas J Silhavy. The essential tension: opposed reactions in bacterial two-component regulatory systems. Trends in microbiology, 1(8):306–310, 1993.
[5] Weihong Hsing, Frank D Russo, Karen K Bernd, and Thomas J Silhavy. Mutations that alter the kinase and phosphatase activities of the two-component sensor EnvZ. Journal of bacteriology, 180(17):4538–4546, 1998.
[6] Eric Batchelor and Mark Goulian. Robustness and the cycle of phosphorylation and dephosphorylation in a two-component regulatory system. Proceedings of the National Academy of Sciences, 100(2):691–696, 2003.
[7] Guy Shinar, Ron Milo, María Rodríguez Martínez, and Uri Alon. Input–output robustness in simple bacterial signaling systems. Proceedings of the National Academy of Sciences, 104(50):19931–19935, 2007.
[8] Uri Alon. An introduction to systems biology: design principles of biological circuits. CRC press, 2019.
[9] Yuval Hart and Uri Alon. The utility of paradoxical components in biological circuits. Molecular cell, 49(2):213–221, 2013.
[10] Jinsu Kim and German Enciso. Absolutely robust controllers for chemical reaction networks. Journal of the Royal Society Interface, 17(166):20200031, 2020.
[11] Fangzhou Xiao and John C Doyle. Robust perfect adaptation in biomolecular reaction networks. In 2018 IEEE conference on decision and control (CDC), pages 4345–4352. IEEE, 2018.
[12] Daniele Cappelletti, Ankit Gupta, and Mustafa Khammash. A hidden integral structure endows absolute concentration robust systems with resilience to dynamical concentration disturbances. Journal of the Royal Society Interface, 17(171):20200437, 2020.
[13] Mustafa H Khammash. Perfect adaptation in biology. Cell Systems, 12(6):509–521, 2021.
[14] Eduardo D Sontag. Adaptation and regulation with signal detection implies internal model. Systems & control letters, 50(2):119–126, 2003.
[15] Ankit Gupta and Mustafa Khammash. Universal structural requirements for maximal robust perfect adaptation in biomolecular networks. Proceedings of the National Academy of Sciences, 119(43):e2207802119, 2022.
[16] Badal Joshi and Anne Shiu. A survey of methods for deciding whether a reaction network is multistationary. “Chemical Dynamics” – special issue of Mathematical Modelling of Natural Phenomena, 10(5):47–67, 2015.
[17] Polly Y Yu and Gheorghe Craciun. Mathematical analysis of chemical reaction systems. Israel Journal of Chemistry, 58(6-7):733–741, 2018.
[18] David G. Schaeffer and John W. Cain. Ordinary differential equations: Basics and beyond. Springer, 2018.
[19] John Guckenheimer and Philip Holmes. Nonlinear oscillations, dynamical systems, and bifurcations of vector fields, volume 42. Springer Science & Business Media, 2013.
[20] Nicolette Meshkat, Anne Shiu, and Angelica Torres. Absolute concentration robustness in networks with low-dimensional stoichiometric subspace. Vietnam Journal of Mathematics, 50:623–651, 2022.
[21] Badal Joshi and Tung D Nguyen. Bifunctional enzyme provides absolute concentration robustness in multisite covalent modification networks. arXiv preprint arXiv:2304.02679, 2023.

Power-engine-load form for dynamic absolute concentration robustness

Abstract

1 Introduction

2 Reaction networks

2.1 Mass action systems

2.2 Dynamic absolute concentration robustness

Definition 2.1.

Definition 2.2.

3 Power-engine-load form of dynamic ACR

3.1 Examples of power-engine-load form in mass action systems

Example 3.1.

Example 3.2.

Example 3.3.

3.2 Convergence results for power-engine-load form dynamic ACR

Theorem 3.4 (Zero load).

Proof.

Example 3.5.

Corollary 3.6.

Proof.

Example 3.7.

Theorem 3.8.

Proof.

4 Applications of power-engine-load form

4.1 Motifs of static and dynamic ACR

4.2 Positive slope motif

Theorem 4.1.

Proof.

Corollary 4.2.

4.3 Negative slope motif

Theorem 4.3.

Proof.

Remark 4.4.

Corollary 4.5.

Proof.

Example 4.6.

Example 4.7.

Remark 4.8.

Theorem 4.9.

Proof.

4.4 Negative slope motif with inflows and outflows

Theorem 4.10.

Proof.

Remark 4.11.

Corollary 4.12.

Proof.

Example 4.13.

Example 4.14.

4.5 A representative of the negative slope motif with stronger robustness properties

Theorem 4.15.

Proof.

Theorem 4.16.

Proof.

Remark 4.17.

Remark 4.18.

Example 4.19.

4.6 Zero slope motif

Theorem 4.20.

Proof.

Remark 4.21.

4.7 Enzyme catalysis with bifunctional enzyme

Theorem 4.22.

Proof.

Theorem 4.23.

Proof.

5 Discussion

Acknowledgments

References

Power-engine-load form for
dynamic absolute concentration robustness