Reaction Network Motifs for Static and 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

Networks with absolute concentration robustness (ACR) have the property that a translation of a coordinate hyperplane either contains all steady states (static ACR) or attracts all trajectories (dynamic ACR). The implication for the underlying biological system is robustness in the concentration of one of the species independent of the initial conditions as well as independent of the concentration of all other species. Identifying network conditions for dynamic ACR is a challenging problem. We lay the groundwork in this paper by studying small reaction networks, those with 2 reactions and 2 species. We give a complete classification by ACR properties of these minimal reaction networks. The dynamics is rich even within this simple setting. Insights obtained from this work will help illuminate the properties of more complex networks with dynamic ACR.
Keywords: Reaction networks, Absolute Concentration Robustness, dynamic ACR, robustness

1 Introduction

Absolute Concentration Robustness (ACR) was introduced in [22] as the mathematical property that every positive steady state of a reaction system coincides in one of the coordinates. We call this property static ACR. Assuming convergence to one of these positive steady states, static ACR ensures that at steady state, the value of one of the measured variables will be independent of the initial value. However, such convergence is not guaranteed even for the simplest systems. The long-term behavior of a system can lead to the system converging to the boundary, or to another attractor such as a limit cycle or not converging at all but diverging to infinity. We introduced dynamic ACR in [13] to account for the global dynamics, and not merely the location of steady states.

A network condition for static ACR is found in [22]: if the network deficiency is one and two non-terminal complexes differ in a single species, then the concentration of that species will have static ACR for all mass action reaction rate constants. Since analyzing global dynamics of a dynamical system is an enormously more complex task than determining the locations of steady states, obtaining network conditions for dynamic ACR is challenging. We make some headway in this direction by studying small reaction networks and organizing our ideas carefully in this setting. In this paper, we study networks with two reactions and two species and make fine distinctions between global dynamical properties related to convergence to an ACR hyperplane. This lays the groundwork for future study of network conditions for larger, more biochemically realistic networks.

The definition of dynamic ACR introduced in [13] is generalized here in two ways. The first is by including a basin of attraction for the ACR hyperplane, which may be the entire positive orthant or a proper subset of it. The second way is by considering a weaker version of attracting hyperplane: trajectories may not converge to the ACR hyperplane, but may only move in the direction of that hyperplane. For static ACR, we introduce a stronger form which requires at least one steady state in each compatibility class that intersects the ACR hyperplane. While each ACR property (static, strong static, dynamic, weak dynamic) is worthy of study by itself, there are connections between them which give a more complete picture. Weak dynamic ACR is a necessary condition for dynamic ACR; similarly, when a positive steady state exists, static ACR is necessary for dynamic ACR (see Theorem 2.12 and Figure 4). Therefore networks with static ACR and weak dynamic ACR serve as candidates for networks with dynamic ACR. We show that all motifs of static ACR and weak dynamic ACR for networks with two reactions and two species can be completely characterized – there are exactly 8 network motifs with static ACR and there are exactly 17 network motifs with weak dynamic ACR and these are depicted in Figure 1 (see Theorems 3.5, 3.6 and 4.1; also see [18] for a network characterization of static ACR).

Each network motif is related to an infinite family of networks, see Section 3.4 for an explanation of how a network maps to a motif. The motifs in Figure 1 (or rather the networks that map to the motifs) have a rich diversity of dynamical properties. Key results in this paper, appearing in Section 3 onwards, deal with identifying such properties. The results are summarized in Figures 13 and 14 – fully annotated versions of Figure 1. The reader is encouraged to begin with at least a glance at the summary theorem (Theorem 5.1) and Figures 13 and 14 to orient themselves towards the objectives of the paper. Specifically, we identify which of the 8 network motifs with static ACR and which of the 17 network motifs with weak dynamic ACR are also dynamic ACR.

(a) The 8 network motifs with static ACR

(b) The 17 network motifs with weak dynamic ACR

Figure 1: Network motifs with output robustness: (Left) All possible two-reaction two-species network motifs with static ACR. (Right) All possible two-reaction two-species network motifs for which the hyperplane

\{x=x^{*}\}

is weakly stable (invariant and weakly attracting) (see Definition 2.9.) The motifs with one-dimensional stoichiometric space are on the horizontal band. The motifs on the circumference require at least two species while the motif in the center can be realized with one or two species.

Infinitely many reaction networks map onto a single network motif – an object that captures two crucial network properties common to the family of networks: (i) orientation of the reactant polytope, (ii) orientation of the reaction vectors relative to the coordinate axes and relative to one another. These two network properties are sufficient for characterizing a wide range of dynamic ACR properties. To understand these network properties, we embed a reaction network in Euclidean space by identifying the stoichiometric coefficients of each species in a complex with a point in Euclidean space. For instance, the reactant and product complex of the reaction $X+Y\to 2Z$ are identified with the points $(1,1,0)$ and $(0,0,2)$ in Euclidean space. The reactant polytope of a reaction network is the convex hull of the reactant complexes in Euclidean space. In the case of a reaction network with only two reactions, a reactant polytope is merely a line segment or a single point. Reactions are fixed “arrows” in Euclidean space. For the case of two species, a reaction may be aligned with the coordinate axes or may be in the interior of a quadrant. Moreover when there are two reactions, multiple relative orientations of the reaction arrows are possible. A network motif is depicted by drawing the orientation of the reactant polytope and the relative orientations of the reactions. A network motif is deeply connected with the dynamics of the reaction system. For one-dimensional systems such as the one in Figure 2(a), the reaction vectors must be parallel to one another and moreover the trajectories in phase space, when superimposed on the reaction network embedding, are simply parallel to the reaction vectors. For two and higher dimensional systems, the relationship is more complicated since the vector field is a linear combination of the reaction vectors with variable coefficients that depend on the position in phase space. The reactant polytope is extremely significant for ACR. All the various forms of ACR properties with at most two species require that the reactant polytope be parallel to one of the coordinate axes, i.e. the polytope should be a horizontal or a vertical line segment. When the reactant polytope is horizontal, the ACR hyperplane is vertical and vice versa. Among static ACR networks in two species, we show that dynamic ACR requires that the reactions points “inwards”, along the reactant polytope. For the weak dynamic ACR networks in two species, we show that dynamic ACR requires that on average, reactions point away from the coordinate axis in the direction perpendicular to the reactant polytope (“up” when the reactant polytope is horizontal). Of course, the above statements are merely suggestive, see Theorem 5.1 for precise statements.

Related work: The idea of minimal reaction networks with a certain dynamical property has a long history. While this work is the first to identify minimal motifs of dynamic ACR, network motifs for other dynamical properties have been studied extensively. For example, networks with two species which have limit cycles have been identified in [21], while [20] studies simple networks with phase transitions giving candidates for multistationarity. In recent years, several classes of minimal multistationary networks called atoms of multistationarity have been identified in [14, 12, 15, 9], and those of oscillations in [5]. Network conditions for static ACR and minimal networks with static ACR have also appeared in recent work [23, 19, 18]. Biochemical implications of static ACR have been studied [17, 11, 10] as well as the implications from a control theory perspective [8]. Stochastic (continuous-time Markov chain) models of reaction networks with the ACR property were studied in [4, 3].

This article is organized as follows. Section 2 introduces different forms of static and dynamic ACR and the relations between them. Section 3 gives a classification of small reaction networks, those containing at most two reactions and two species; here the focus is on mass-conserving reaction networks. Section 4 continues the classification, extended to networks where the trajectories may go to infinity or to the boundary. Section 5 summarizes the main classification results with the statement of a larger theorem that encapsulates numerous results from the previous sections.

2 Forms of ACR

Background and terminology

We use standard notation and terminology for reaction networks and mass action systems. A quick summary is given here, see for instance [15, 26] for further details. Upper case letters ( $X,Y,Z,A,B$ ) are used for species participating in reactions and the corresponding lower case letters ( $x,y,z,a,b$ ) for their concentrations which are time-varying quantities. An example of a reaction is $X+Y\to 2Z$ , where $X+Y$ is referred to as the source complex, while $2Z$ is the product complex. The rate of any given reaction is a nonnegative-valued function of species concentrations. In the case of mass action kinetics, the rate is proportional to the product of reactant concentrations taken with multiplicity. The proportionality constant, called the reaction rate constant, is placed adjacent to the reaction arrow, as follows: $X+Y\xrightarrow{k}2Z$ . The rate of this reaction under mass action kinetics is $kxy$ . The reaction vector for this reaction is the difference between the product complex and the source complex, i.e. $2Z-(X+Y)$ , which under a choice of a standard coordinate basis is written as $(-1,-1,2)$ . A reaction network is a nonempty set of reactions, such that every species participates in at least one reaction, and none of the reaction vectors is the zero vector. The stoichiometric subspace of a reaction network is the subspace spanned by the set of reaction vectors of the reaction network.

We use $\mathcal{G}$ to denote a reaction network and $K$ to denote a specific choice of mass action kinetics for $\mathcal{G}$ , so that $(\mathcal{G},K)$ is a mass action dynamical system. Throughout the paper, we depict the Euclidean embedding of a reaction network next to the dynamics in the phase plane. Figure 2(a) is an example of a Euclidean embedding: the reactions $B\to A$ and $A+B\to 2B$ are shown as red arrows, each arrow originating at the reactant complex and terminating at the product complex. We also depict the reactant polytope, defined as the convex hull of the reactant complexes. In the example in Figure 2(a), the set of reactant complexes is $\{B,A+B\}$ and so the reactant polytope, shown in green, is a line segment joining the two complexes.

Throughout this paper, we consider a dynamical system $\mathcal{D}$ defined by $\dot{x}=f(x)$ with $x\in\mathbb{R}^{n}_{\geq 0}$ for which $\mathbb{R}^{n}_{\geq 0}$ is forward invariant. $x\in\mathbb{R}^{n}_{\geq 0}$ is a steady state of $\mathcal{D}$ if $f(x)=0$ .

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

The notation $\mathcal{H}[i,a_{i}^{*}]\coloneqq\{x\in\mathbb{R}^{n}_{>0}:x_{i}=a_{i}^{*}\}$ is reserved for the hyperplane parallel to a coordinate hyperplane and restricted to the positive orthant. When variables are labelled without indices, for instance as $x,y,z$ etc., we use the notation $\mathcal{H}[x,x^{*}]\coloneqq\{x\in\mathbb{R}^{n}_{>0}:x=x^{*}\}$ . Even though the two notations are slightly inconsistent, there is no possibility of confusion.

Static ACR and dynamic ACR in a real dynamical system were defined in [13]. We repeat these definitions here:

Definition 2.2.

•

$\mathcal{D}$ is a static ACR system if $\mathcal{D}$ has a positive steady state and there is an $i\in\{1,\ldots,n\}$ and a positive $a_{i}^{*}\in\mathbb{R}_{>0}$ such that any positive steady state $x\in\mathbb{R}^{n}_{>0}$ is in the hyperplane $\{x_{i}=a_{i}^{*}\}$ . Any such $x_{i}$ and $a_{i}^{*}$ is a static ACR variable and its static ACR value, respectively. $\mathcal{H}[i,a_{i}^{*}]$ is the static ACR hyperplane.
•

$\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 $x(0)\in\mathbb{R}^{n}_{>0}$ that is compatible with $\{x\in\mathbb{R}^{n}_{>0}~|~x_{i}=a_{i}^{*}\}$ , a unique solution to $\dot{x}=f(x)$ exists up to some maximal $T_{0}(x(0))\in(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. $\mathcal{H}[i,a_{i}^{*}]$ is the dynamic ACR hyperplane.

Now we define strong and weak forms of each type of ACR. All forms: static ACR, strong static ACR, dynamic ACR and weak dynamic ACR are marked by the presence of a privileged hyperplane parallel to a coordinate hyperplane $\{x\in\mathbb{R}^{n}_{>0}:x_{i}=a_{i}^{*}\}$ called the ACR hyperplane. This hyperplane either contains all the steady states (in the static case) or is the unique attractor for all relevant trajectories (in the dynamic case). Moreover, each form of ACR has a set $\Omega$ associated with it called the basin of ACR.

2.1 Static forms of ACR

We generalize the definition of static ACR from [13] to allow for arbitrary basin sets $\Omega$ .

Definition 2.3.

The variable $x_{i}$ , where $i\in\{1,\ldots,n\}$ , has static ACR w.r.t. $\Omega\subseteq\mathbb{R}_{\geq 0}^{n}$ if there is an $a_{i}^{*}>0$ such that the following hold:

1.

$f(x)=0$ for some positive $x\in\Omega$ ,
2.

for any $x\in\Omega$ such that $f(x)=0$ , $x_{i}=a_{i}^{*}$ .

In this case, the static ACR value of $x_{i}$ is $a_{i}^{*}$ and $\mathcal{H}[i,a_{i}^{*}]\coloneqq\{x\in\mathbb{R}^{n}_{>0}:x_{i}=a_{i}^{*}\}$ is the static ACR hyperplane.

Remark 2.4.

When $\Omega=\mathbb{R}^{n}_{>0}$ , we simply say that a variable has “static ACR” instead of “static ACR w.r.t. $\mathbb{R}^{n}_{>0}$ ”. Note that boundary steady states are conventionally excluded from consideration, in other words, we do not require that $\Omega=\mathbb{R}^{n}_{\geq 0}$ for static ACR, merely that $\Omega=\mathbb{R}^{n}_{>0}$ .

By strong static ACR, we mean the property that every compatibility class that intersects the static ACR hyperplane contains at least one positive steady state. Most commonly studied motifs and biochemical systems with static ACR are strong static ACR as well.

Definition 2.5.

The variable $x_{i}$ , where $i\in\{1,\ldots,n\}$ , has strong static ACR w.r.t. $\Omega\subseteq\mathbb{R}_{\geq 0}^{n}$ if there is an $a_{i}^{*}>0$ such that the following hold:

1.

$x_{i}$ has static ACR w.r.t. $\Omega\subseteq\mathbb{R}_{\geq 0}^{n}$ with value $a_{i}^{*}$ ,
2.

for any $y\in\mathcal{H}[i,a_{i}^{*}]\cap\Omega$ there is a $z\in\Omega\cap\mathbb{R}^{n}_{>0}$ such that $y-z\in\SS_{f}$ and $f(z)=0$ .

In this case, the strong static ACR value of $x_{i}$ is $a_{i}^{*}$ and $\mathcal{H}[i,a_{i}^{*}]$ is the strong static ACR hyperplane.

The concentration of $A$ in the network in Figure 2 has strong static ACR w.r.t. $\mathbb{R}^{2}_{>0}$ because every point on the hyperplane $\mathcal{H}[a,a^{*}]=\{(a,b)\in\mathbb{R}^{2}_{>0}:a=a^{*}\}$ (the green vertical line) is a steady state. So any compatibility class with $a(0)+b(0)>a^{*}$ intersects the hyperplane $\mathcal{H}[a,a^{*}]$ and therefore has a positive steady state. Every network motif with static ACR studied in this paper also has strong static ACR, as we will show in Theorem 3.5.

Refer to caption — (a) Reaction network embedded in Euclidean plane

2.2 Dynamic forms of ACR

We generalize the definition of dynamic ACR from [13] in two ways: (i) we allow for arbitrary basin sets $\Omega$ , and (ii) we define a weaker form for which the ACR hyperplane is weakly attracting.

Definition 2.6.

The variable $x_{i}$ , where $i\in\{1,\ldots,n\}$ , has dynamic ACR w.r.t. $\Omega\subseteq\mathbb{R}_{\geq 0}^{n}$ if there is an $a_{i}^{*}>0$ such that for any initial value $x(0)$ in $\Omega$ , a unique solution to $\dot{x}=f(x)$ exists up to some maximal $T_{0}(x(0))\in(0,\infty]$ , and $\lim_{t\to T_{0}}x_{i}(t)=a_{i}^{*}$ . We say that $x_{i}$ has dynamic ACR value $a_{i}^{*}$ . Moreover, we say that the ACR hyperplane $\mathcal{H}[i,a_{i}^{*}]$ is an attractor for $\Omega$ and that $\Omega$ is a basin of attraction of $\mathcal{H}[i,a_{i}^{*}]$ .

Weak dynamic ACR is the notion that the ACR variable converges to a value that is not further from the ACR value than the initial distance. In other words, all initial conditions are attracted towards the ACR hyperplane even if they fail to reach there.

Definition 2.7.

The variable $x_{i}$ , where $i\in\{1,\ldots,n\}$ , has weak dynamic ACR w.r.t. $\Omega\subseteq\mathbb{R}^{n}_{\geq 0}$ if there is an $a_{i}^{*}>0$ such that for any initial value $x(0)$ in $\Omega$ , a unique solution to $\dot{x}=f(x)$ exists up to some maximal $T_{0}(x(0))\in(0,\infty]$ , $\lim_{t\to T_{0}}x_{i}(t)$ exists and $\lim_{t\to T_{0}}\absolutevalue{x_{i}(t)-a_{i}^{*}}\leq\absolutevalue{x_{i}(0)-a_{i}^{*}}$ , with strict inequality when $x_{i}(0)\neq a_{i}^{*}$ . We say that $x_{i}$ has weak dynamic ACR value $a_{i}^{*}$ . Moreover, we say that the ACR hyperplane $\mathcal{H}[i,a_{i}^{*}]$ is a weak attractor for $\Omega$ and $\Omega$ is a weak basin of attraction of $\mathcal{H}[i,a_{i}^{*}]$ .

Remark 2.8.

We do not require the (weak) ACR hyperplane is invariant, merely that it is a (weak) attractor. A weak attractor may fail to be an attractor because trajectories reach the boundary or diverge to infinity before reaching the ACR hyperplane.

In analogy with the definition of stability of a steady state, we define the following.

Definition 2.9.

•

The hyperplane $\mathcal{H}[i,a_{i}^{*}]$ is stable w.r.t. $\Omega$ if $\mathcal{H}[i,a_{i}^{*}]$ is an attractor for $\Omega$ and all trajectories with initial value in $\Omega$ move monotonically towards $\mathcal{H}[i,a_{i}^{*}]$ .
•

The hyperplane $\mathcal{H}[i,a_{i}^{*}]$ is weakly stable w.r.t. $\Omega$ if $\mathcal{H}[i,a_{i}^{*}]$ is a weak attractor for $\Omega$ and all trajectories with initial value in $\Omega$ move monotonically towards $\mathcal{H}[i,a_{i}^{*}]$ .

Remark 2.10.

Stability (weak or not) of $\mathcal{H}[i,a_{i}^{*}]$ implies that $\mathcal{H}[i,a_{i}^{*}]$ is invariant.

Example 2.11.

(Weakly stable but not stable hyperplane) Consider the mass action system of the following reaction network (see also Figure 3(a)):

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

The mass action system of ODEs is:

\displaystyle\dot{a}=b(k_{2}-2k_{1}a^{2}),\quad\dot{b}=-b(k_{2}-k_{1}a^{2}).

(2.1)

It is clear from $\dot{a}$ that for the duration of time that $b(t)$ is positive, $a(t)$ approaches the value $\sqrt{k_{2}/(2k_{1})}$ . This implies that $a$ is a weak dynamic ACR variable with a weak ACR value of $\sqrt{k_{2}/(2k_{1})}$ . Note however that $\dot{a}+\dot{b}=-k_{1}a^{2}b$ is negative everywhere in the positive orthant, and so we expect $b(t)$ to converge to $0$ for most initial values. In fact, this is the case for every initial value that is not on the weak ACR hyperplane, as a consequence of Theorem 4.2.

2.3 Relations between different forms of ACR

It is clear from the definitions that if $x_{i}$ has dynamic ACR w.r.t. $\Omega$ with value $a_{i}^{*}$ , then $x_{i}$ has weak dynamic ACR w.r.t. $\Omega$ with the same value $a_{i}^{*}$ . Similarly strong static ACR w.r.t. $\Omega$ implies static ACR w.r.t. $\Omega$ . The static forms of ACR are also related to the dynamic forms under some mild assumptions of existence of steady states.

Theorem 2.12.

Suppose that $x_{i}$ is weak dynamic ACR w.r.t. $\Omega$ with value $a_{i}^{*}$ .

1.

Suppose there is an $x\in\Omega$ such that $f(x)=0$ . Then $x_{i}$ is static ACR w.r.t. $\Omega$ .
2.

Suppose for any $y\in\mathcal{H}[i,a_{i}^{*}]\cap\Omega$ there is a positive $x\in\Omega$ such that $y-x\in\SS_{f}$ and $f(x)=0$ . Then $x_{i}$ is strong static ACR w.r.t. $\Omega$ .

Proof.

For both cases, the static ACR property follows from observing that there cannot be any positive steady state in $\Omega\setminus\mathcal{H}[i,a_{i}^{*}]$ because such a steady state violates the weak dynamic ACR hypothesis. The additional implication of strong static ACR in the second case is immediate from the definition. ∎

These relations are portrayed in Figure 4.

Figure 4: (Relations between static ACR and dynamic ACR) The implications in black follow from the definitions. The implication in cyan requires the additional assumption of existence of a positive steady state. The implication in magenta requires the additional assumption of existence of a positive steady state within each compatibility class which has a nonempty intersection with the ACR hyperplane

[i,a_{i}^{*}]

2.4 Some basins of interest for all forms of ACR

We discuss some basins of natural interest and the relations between them. A basin of ACR applies to any of the forms of ACR (static, strong static, dynamic, weak dynamic) discussed earlier.

Definition 2.13.

Let $\mathcal{P}\in\{$ static, strong static, dynamic, weak dynamic $\}$ . We define various basin types as follows.

(i)

Full basin $\mathcal{P}$ -ACR occurs when $\Omega=\mathbb{R}^{n}_{>0}$ . Full basin static ACR is simply referred to as static ACR.
(ii)

Subspace $\mathcal{P}$ -ACR occurs when $\Omega=(\mathcal{H}[i,a_{i}^{*}]+\SS)\cap\mathbb{R}^{n}_{>0}$ for some subspace $\SS$ of $\SS_{f}$ such that $\SS\not\subseteq e_{i}^{\perp}$ . We say full space $\mathcal{P}$ -ACR if we have subspace $\mathcal{P}$ -ACR with $\SS=\SS_{f}$ . Full space dynamic ACR is simply referred to as dynamic ACR.
(iii)

Neighborhood $\mathcal{P}$ -ACR occurs when $\Omega$ is a neighborhood $\mathcal{H}[i,a_{i}^{*}]$ .

Suppose there are some $M_{j}>0$ for all $j\neq i$ such that the neighborhood of the set $\{x_{i}=a_{i}^{*},x_{j}>M_{j}:j\neq i\}$ is a basin of $\mathcal{P}$ -ACR. Then we say that $x_{i}$ has almost neighborhood $\mathcal{P}$ -ACR.
(iv)

Cylinder $\mathcal{P}$ -ACR occurs when $\Omega$ is a cylinder of $\mathcal{H}[i,a_{i}^{*}]$ . A cylinder of $\mathcal{H}[i,a_{i}^{*}]$ is the set of points $\{\absolutevalue{x_{i}-a_{i}^{*}}<\delta^{*}\}$ for some $\delta^{*}>0$ .

We define almost cylinder $\mathcal{P}$ -ACR when a basin of $\mathcal{P}$ -ACR is a cylinder of some set $\{x_{i}=a_{i}^{*},x_{j}>M_{j}:j\neq i\}$ for some $M_{j}>0$ for all $j\neq i$ .
(v)

Null $\mathcal{P}$ -ACR occurs when $\Omega=\mathcal{H}[i,a_{i}^{*}]$ . Non-null $\mathcal{P}$ -ACR occurs when $\Omega\setminus\mathcal{H}[i,a_{i}^{*}]\neq\varnothing$ .

Figure 5: (Relations between basin types.) Let

\mathcal{P}\in\{

static, strong static, dynamic, weak dynamic

\}

. The basin type implications are based on the observation that if

\Omega\subseteq\Omega^{\prime}

then

\mathcal{P}

-ACR w.r.t.

\Omega^{\prime}

implies

\mathcal{P}

-ACR w.r.t.

\Omega

Neighborhood ACR, cylinder ACR, almost neighborhood ACR and almost cylinder ACR are local forms of ACR. Full basin and subspace ACR are non-local forms of ACR. Null $\mathcal{P}$ -ACR is considered neither local nor non-local. Conversely, all local and non-local forms of $\mathcal{P}$ -ACR are non-null.

It’s clear that if $\Omega\subseteq\Omega^{\prime}$ then $\mathcal{P}$ -ACR w.r.t. $\Omega^{\prime}$ implies $\mathcal{P}$ -ACR w.r.t. $\Omega$ . Certain relations between non-local forms: (full basin $\mathcal{P}$ -ACR $\implies$ full space $\mathcal{P}$ -ACR $\implies$ subspace $\mathcal{P}$ -ACR $\implies$ null $\mathcal{P}$ -ACR) and between local forms: (cylinder $\mathcal{P}$ -ACR $\implies$ neighborhood $\mathcal{P}$ -ACR $\implies$ null $\mathcal{P}$ -ACR & almost neighborhood $\mathcal{P}$ -ACR) and (cylinder $\mathcal{P}$ -ACR $\implies$ almost cylinder $\mathcal{P}$ -ACR $\implies$ almost neighborhood $\mathcal{P}$ -ACR) follow from the definitions. Certain local forms are related to subspace $\mathcal{P}$ -ACR as the following theorem shows.

Theorem 2.14.

Let $\mathcal{D}$ be a dynamical system which is subspace $\mathcal{P}$ -ACR. Then $\mathcal{D}$ is (i) almost cylinder $\mathcal{P}$ -ACR and (ii) neighborhood $\mathcal{P}$ -ACR.

Proof.

Let $\SS$ be a subspace of the stoichiometric space $\SS_{f}$ such that $x_{1}$ is a subspace $\mathcal{P}$ -ACR variable with value $a_{1}^{*}$ . Let $v\in\SS\cap\left(\mathbb{R}^{n}\setminus e_{1}^{\perp}\right)$ , i.e. $v=(v_{1},\ldots,v_{n})$ is some vector with $v_{1}\neq 0$ . Let $\varepsilon\in(0,a_{1}^{*})$ and define the almost cylinder neighborhood of $\mathcal{H}[i,a_{i}^{*}]$ :

\Omega_{\varepsilon}\coloneqq\left\{z\in\mathbb{R}^{n}_{>0}:\absolutevalue{z_{1}-a_{1}^{*}}<\varepsilon,~z_{2}>\varepsilon\frac{\absolutevalue{v_{2}}}{\absolutevalue{v_{1}}},\ldots,z_{n}>\varepsilon\frac{\absolutevalue{v_{n}}}{\absolutevalue{v_{1}}}\right\}.

To see $x_{1}$ is $\mathcal{P}$ -ACR w.r.t. $\Omega_{\varepsilon}$ for every $\varepsilon\in(0,a_{1}^{*})$ , we only need to show that $\Omega_{\varepsilon}$ is contained in $\mathcal{H}[1,a_{1}^{*}]+{\rm span}\{v\}$ which in turn is clearly contained in $\mathcal{H}[1,a_{1}^{*}]+\SS$ . Indeed, let $z\in\Omega_{\varepsilon}$ . Then $z\in\mathcal{H}[1,a_{1}^{*}]+{\rm span}\{v\}$ if there is a $\beta\in\mathbb{R}$ such that $z-\beta v\in\mathcal{H}[1,a_{1}^{*}]$ . Let $\beta\coloneqq(z_{1}-a_{1}^{*})/v_{1}$ . Then

z-\beta v=\left(a_{1}^{*},z_{2}-\frac{z_{1}-a_{1}^{*}}{v_{1}}v_{2},\ldots,z_{n}-\frac{z_{1}-a_{1}^{*}}{v_{1}}v_{n}\right).

For $j\in\{2,\ldots,n\}$ ,

\absolutevalue{\frac{z_{1}-a_{1}^{*}}{v_{1}}v_{j}}=\absolutevalue{z_{1}-a_{1}^{*}}\frac{\absolutevalue{v_{j}}}{\absolutevalue{v_{1}}}<\varepsilon\frac{\absolutevalue{v_{j}}}{\absolutevalue{v_{1}}}<z_{j},

which shows that all components of $z-\beta v$ are positive and the first component is $a_{1}^{*}$ . This proves the almost cylinder property where the almost cylinder is $\Omega_{\varepsilon}$ for any $\varepsilon\in(0,a_{1}^{*})$ .

Let $\Omega\coloneqq\cup_{\varepsilon>0}\Omega_{\varepsilon}$ . It is clear that $\Omega$ is contained in $\mathcal{H}[1,a_{1}^{*}]+\SS$ and that $\Omega$ is a neighborhood of $\mathcal{H}[1,a_{1}^{*}]$ . This shows that $x_{1}$ is neighborhood $\mathcal{P}$ -ACR. ∎

The relations between $\mathcal{P}$ -ACR with different basin types are depicted in Figure 5.

2.5 Some examples of systems with local dynamic ACR

(Almost cylinder ACR/Neighborhood ACR/Weak cylinder ACR) Consider the reaction network shown in Figure 6(a).

(a)

For the mass action system defined by $k_{1}=6,k_{2}=11,k_{3}=6,k_{4}=1$ , the positive steady states form three distinct rays $a=1$ , $a=2$ and $a=3$ . Therefore, $a$ is a cylinder static ACR variable with multiple ACR values $\{1,2,3\}$ . It is easily checked (see simulated trajectories in Figure 6(b)), that only $a=2$ is locally stable within each compatibility class. The maximal cylinder which forms the basin of attraction for $\{a=2\}$ has radius $1$ .
(b)

For the mass action system defined by $k_{1}=1,k_{2}=3,k_{3}=3,k_{4}=1$ , all positive steady states lie on $\{a=1\}$ . Moreover, the positive steady states are repelling. It follows that $a$ is a (global) static ACR variable but not a local dynamic ACR variable.

(c)

For the general case, the system of mass action ODEs is

\displaystyle\dot{a}=-ab(k_{1}-k_{2}a+k_{3}a^{2}-k_{4}a^{3}),\quad\dot{b}=ab(k_{1}-k_{2}a+k_{3}a^{2}-k_{4}a^{3}).

The univariate polynomial $k_{1}-k_{2}a+k_{3}a^{2}-k_{4}a^{3}$ must have at least one positive, real zero. Moreover, when there are multiple positive zeros, say $a^{*}$ and $a^{**}$ , clearly there are non-intersecting cylinders that contain the sets $\{a=a^{*}\}$ and $\{a=a^{**}\}$ . This shows that $A$ is a local static ACR species, i.e. for any choice of mass action kinetics, $a$ is a local static ACR variable. In some cases, $a$ may be a global static ACR variable. However, $A$ is not a local dynamic ACR species as there is a choice of rate constants for which $a$ is not a local dynamic ACR variable. In other words, the reaction network has capacity for local dynamic ACR and is local static ACR (see Definition 3.1 for meaning of ‘capacity’).

2.

(Multiple ACR values in almost cylinder ACR) Consider the reaction network shown in Figure 6(a).

(a) Reaction network embedded in Euclidean plane

(b) Trajectories in phase plane

Figure 7: The mass action system of the reaction network ( $A+B\xrightarrow{k_{2}}2B,~B\xrightarrow{k_{1}}A,~3A+B\xrightarrow{k_{4}}2A+2B,~2A+B\xrightarrow{k_{3}}3A$ ). For the rate constants $k_{1}=6,k_{2}=11,k_{3}=6,k_{4}=1$ , $a$ is a local static and dynamic ACR variable. There are three local static ACR values, $a=1,a=2$ , and $a=3$ , of which $a=1$ and $a=3$ are local dynamic ACR values.

(Neighborhood ACR but not almost cylinder ACR) Consider the mass action system shown below:

	$\displaystyle A+B\xrightarrow{1}2A,$	$\displaystyle\quad 2A+B\stackrel{{\scriptstyle[}}{{1}}]{1}{\rightleftarrows}A+2B,$
	$\displaystyle 2A+2B\xrightarrow{2}A+3B,$	$\displaystyle\quad 3A+2B\xrightarrow{1}4A+B.$		(2.2)

The resulting mass action system when taken with inflows is

	$\displaystyle\dot{a}$	$\displaystyle=ab(1-a)(1-(a-1)b)+g_{a},$
	$\displaystyle\dot{b}$	$\displaystyle=-ab(1-a)(1-(a-1)b)+g_{b}.$

Some trajectories are shown in Figure 8 for the case of $g_{a}=0.1$ and $g_{b}=0$ . Any initial value to the left of the ACR hyperplane converges to the ACR hyperplane, but for every initial $a$ value to the right of the ACR hyperplane, there is a large enough $b$ value such that the trajectory does not converge to the ACR value.

\suspend

enumerate

2.6 An example of a system with subspace (but not full space) dynamic ACR

The following example illustrates the need for defining subspace ACR.

\resume

enumerate

4.

(Subspace ACR and Cylinder ACR) Consider the following mass action system:

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

which defines the ODE system:

$\displaystyle\dot{a}=b(k_{2}-k_{1}a),\quad\dot{b}=-b(k_{2}-2k_{1}a).$

As long as $b(t)$ remains positive, $a(t)$ will move towards $a^{*}=k_{2}/k_{1}$ . Since the stoichiometric subspace is all of $\mathbb{R}^{2}$ , dynamic ACR requires every positive initial value to converge to $a^{*}$ . This condition is not satisfied since, as trajectories in bottom left corner of Figure 9(b) show, there exist initial conditions which converge to the $b=0$ boundary and not to the ACR hyperplane. However, if we define $\Omega$ to be $\{a=a^{*}\}+{\rm span}\{(-1,1)\}$ , then $a$ is dynamic ACR w.r.t. $\Omega$ . Thus, $a$ is subspace (dynamic) ACR.

(a) Reaction network embedded in Euclidean plane

(b) Trajectories in phase plane

Figure 9: The mass action system $\{A+B\xrightarrow{k_{1}}3B,B\xrightarrow{k_{2}}A\}$ is subspace dynamic ACR. For any choice of rate constants, all initial values that are $\SS$ -compatible, where $\SS={\rm span}\{(-1,1)\}$ with the ACR hyperplane converge to the ACR hyperplane.

3 Classification of Minimal Static and Dynamic ACR Networks

We consider networks of small size, ones with at most two reactions and at most two species. For such networks, we can catalogue many ACR properties. Some of these networks have archetypal ACR dynamics. The study of minimal, archetypal motifs is valuable because it may reveal the underlying principles at play in the dynamics of more complex networks.

Definition 3.1.

Suppose that $(\mathcal{G},K)$ is a mass action system resulting from the reaction network $\mathcal{G}$ , where $K$ denotes the specific choice of mass action rate constants. Let $\mathcal{P}\in\{$ static, strong static, dynamic, weak dynamic, wide basin dynamic, narrow basin dynamic, full basin dynamic $\}$ .

•

We say that $\mathcal{G}$ has capacity for $\mathcal{P}$ -ACR if there is a $K$ such that the mass action system $(\mathcal{G},K)$ is a $\mathcal{P}$ -ACR system.
•

We say that $\mathcal{G}$ is a $\mathcal{P}$ -ACR network if $(\mathcal{G},K)$ is a $\mathcal{P}$ -ACR system for all choices of $K$ .
•

We say a species $X$ in a network $\mathcal{G}$ is a $\mathcal{P}$ -ACR species if the concentration of $X$ is a $\mathcal{P}$ -ACR variable in $(\mathcal{G},K)$ for all choices of $K$ .

3.1 Static and dynamic ACR for reaction networks with only 1 reaction or only 1 species.

Theorem 3.2 (Static and Dynamic ACR in one-reaction networks).

A network with $n\geq 1$ species and only one reaction is neither static nor dynamic ACR for any choice of rate constants.

Proof.

A network with only one reaction has no positive steady states and is therefore not static ACR. Such a network is also not dynamic ACR since for each $i\in\{1,\ldots,n\}$ , $\dot{x}_{i}$ is either strictly positive in the entire positive orthant, or strictly negative in the entire positive orthant, or identically zero in the entire positive orthant. Therefore, either $x_{i}$ goes to infinity, to zero, or $\dot{x}_{i}\equiv 0$ . In every case $x_{i}$ fails to be a dynamic ACR variable. ∎

Theorem 3.3 (Static ACR in one-species networks).

Let $\mathcal{G}$ be a network with $1$ species and arbitrary number of reactions. The following are equivalent.

A1.

$\mathcal{G}$ has the capacity for static ACR.
A2.

$(\mathcal{G},K)$ has a unique positive steady state for some $K$ .

The following are equivalent.

B1.

$\mathcal{G}$ is static ACR.
B2.

$(\mathcal{G},K)$ has a unique positive steady state for every choice of $K$ .

Proof.

The results follow immediately from basic properties of one dimensional dynamical systems. ∎

Theorem 3.4 (Dynamic ACR in one-species networks).

Let $\mathcal{G}$ be a network with $1$ species and arbitrary number of reactions. The following are equivalent.

A1.

$\mathcal{G}$ has the capacity for weak dynamic ACR.
A2.

$\mathcal{G}$ has the capacity for dynamic ACR.
A3.

$(\mathcal{G},K)$ has a unique positive steady state for some $K$ , and this steady state is globally attracting.
A4.

$(\mathcal{G},K)$ is full basin dynamic ACR for some $K$ .

The following are equivalent.

B1.

$\mathcal{G}$ is weak dynamic ACR.
B2.

$\mathcal{G}$ is dynamic ACR.
B3.

For every choice of $K$ , $(\mathcal{G},K)$ has a unique positive steady state, and this steady state is globally attracting.
B4.

$(\mathcal{G},K)$ is full basin dynamic ACR for every $K$ .

Proof.

The results follow immediately from basic properties of one dimensional dynamical systems. ∎

Many examples of one-species networks along with their ACR properties are shown in Table 1.

Network Arrow Diagram Capacity for static ACR? Is network static ACR? Capacity for dynamic ACR? Is network dynamic ACR? $0\to A$ $\longrightarrow$ No No No No $0\rightleftarrows A$ $\longrightarrow,\longleftarrow$ Yes Yes Yes Yes $0\to A,2A\to 3A$ $\longrightarrow,\longrightarrow$ No No No No $0\to A,2A\leftarrow 3A$ $\longrightarrow,\longleftarrow$ Yes Yes Yes Yes $0\leftarrow A,2A\to 3A$ $\longleftarrow,\longrightarrow$ Yes Yes No No $0\leftarrow A,2A\leftarrow 3A$ $\longleftarrow,\longleftarrow$ No No No No $2A\rightleftarrows 3A$ $\longrightarrow,\longleftarrow$ Yes Yes Yes Yes $0\rightleftarrows A,2A\to 3A$ $\longrightarrow,\longleftarrow,\longrightarrow$ Yes No No No $0\rightleftarrows A,2A\leftarrow 3A$ $\longrightarrow,\longleftarrow,\longleftarrow$ Yes Yes Yes Yes $0\to A,2A\rightleftarrows 3A$ $\longrightarrow,\longrightarrow,\longleftarrow$ Yes Yes Yes Yes $0\leftarrow A,2A\rightleftarrows 3A$ $\longleftarrow,\longrightarrow,\longleftarrow$ Yes No No No $0\rightleftarrows A,~~2A\rightleftarrows 3A$ $\longleftarrow,\longrightarrow,\longleftarrow,\longrightarrow$ Yes No Yes No

Table 1: Subnetworks of

0\rightleftarrows A,~~2A\rightleftarrows 3A

show diverse behaviors when catalogued according to capacities for static and dynamic ACR and according to whether the network is static or dynamic ACR. The diversity illustrates the range of possibilities even for one species networks.

Figure 10: A reaction network with two reactions and at most two species

X

and

Y

can be depicted as a pair of arrows embedded in the Euclidean plane

\mathbb{R}^{2}_{\geq 0}

. The red arrows depict reactions and the green line segment joining the two source complexes is the reactant polytope of a network with two reactions. The arrow from

(a_{1},b_{1})

(\widetilde{a}_{1},\widetilde{b}_{1})

portrays the reaction

a_{1}X+b_{1}Y\to\widetilde{a}_{1}X+\widetilde{b}_{1}Y

. The label

k_{1}

is the mass action rate constant of this reaction. The form of ACR, static or dynamic, as well as basin type can be decided based on the geometry of the three objects appearing in the Figure: the reactant polytope and the reaction arrows.

3.2 Reaction networks with 2 reactions and 2 or fewer species: Notation

We now classify reaction networks with two reactions and two species. We start by defining notation that will be used in the rest of the section. See Figure 10 for a geometric rendering of a reaction network, and how it relates to the notation. Let $\mathcal{G}$ be a reaction network with at most 2 species ( $X,Y$ ) and the following two reactions:

\displaystyle a_{1}X+b_{1}Y\xrightarrow{k_{1}}\widetilde{a}_{1}X+\widetilde{b}_{1}Y,\quad a_{2}X+b_{2}Y\xrightarrow{k_{2}}\widetilde{a}_{2}X+\widetilde{b}_{2}Y,

(3.1)

where $a_{i},b_{i},\widetilde{a}_{i},\widetilde{b}_{i}\in\mathbb{R}_{\geq 0}$ and $(\widetilde{a}_{i},\widetilde{b}_{i})\neq(a_{i},b_{i})$ for $i\in\{1,2\}$ . Although stoichiometric coefficients are usually integers, we allow real values here since the results remain unchanged under this generality. The labels $k_{1}$ and $k_{2}$ are mass action reaction rate constants, and are therefore positive reals. Let

\displaystyle\SS={\rm span}\left\{v_{1}\coloneqq\begin{pmatrix}\widetilde{a}_{1}-a_{1}\\ \widetilde{b}_{1}-b_{1}\end{pmatrix},v_{2}\coloneqq\begin{pmatrix}\widetilde{a}_{2}-a_{2}\\ \widetilde{b}_{2}-b_{2}\end{pmatrix}\right\}

(3.2)

be the stoichiometric subspace of $\mathcal{G}$ . The mass action dynamical system $(\mathcal{G},K)$ explicitly is:

	$\displaystyle\dot{x}$	$\displaystyle=k_{1}(\widetilde{a}_{1}-a_{1})x^{a_{1}}y^{b_{1}}+k_{2}(\widetilde{a}_{2}-a_{2})x^{a_{2}}y^{b_{2}}$
	$\displaystyle\dot{y}$	$\displaystyle=k_{1}(\widetilde{b}_{1}-b_{1})x^{a_{1}}y^{b_{1}}+k_{2}(\widetilde{b}_{2}-b_{2})x^{a_{2}}y^{b_{2}}.$		(3.3)

3.3 Static ACR for reaction networks with 2 reactions and 2 or fewer species.

Theorem 3.5 (Static ACR in networks with 2 reactions and 2 or fewer species).

Let $\mathcal{G}$ be as in (3.1)-(3.2). The following are equivalent:

1.

$\mathcal{G}$ has the capacity for static ACR.
2.

$\mathcal{G}$ is static ACR.
3.

$\mathcal{G}$ is strong static ACR.
4.
- •
  
  the two source complexes are different: $(a_{1},b_{1})\neq(a_{2},b_{2}),$
- •
  
  the source complexes share a common coordinate: $(a_{2}-a_{1})(b_{2}-b_{1})=0$ , and
- •
  
  reaction vectors are negative scalar multiples of each other: $v_{1}=-\mu v_{2}\mbox{ for some }\mu>0,$ (in particular $\dim(\SS)=1$ ).

Furthermore, the following hold when $\mathcal{G}$ has 2 species and is static ACR.

•

Either $X$ or $Y$ , but not both, is a static ACR species.
•

$X$ is an ACR species if $a_{2}\neq a_{1}$ . The variable $x$ has the static ACR value $\displaystyle\left(k_{2}/(\mu k_{1})\right)^{\frac{1}{a_{1}-a_{2}}}$ .
•

$Y$ is an ACR species if $b_{2}\neq b_{1}$ . The variable $y$ has the static ACR value $\displaystyle\left(k_{2}/(\mu k_{1})\right)^{\frac{1}{b_{1}-b_{2}}}$ .

Figure 11: Motifs that do not have the capacity for static ACR. A reaction network does not have the capacity for static ACR if the source complexes of the two reactions are the same, or if both coordinates of the source complexes are different, or if the reaction vectors do not point in opposite directions. See Theorem 3.5 for precise conditions.

Proof.

( ${\it 3}\implies{\it 2}\implies{\it 1}$ ) holds by definition. We now show that ( ${\it 1}\implies{\it 4}$ ). Suppose that $v_{1}\neq-\mu v_{2}$ for any $\mu>0$ , then there are no positive steady states for any choice of mass action rate constants and so $\mathcal{G}$ does not have the capacity for static ACR.

Now assume that $v_{1}=-\mu v_{2}$ for some $\mu>0$ . Then the mass action system is

\displaystyle\dot{x}=(\widetilde{a}_{1}-a_{1})\left(k_{1}x^{a_{1}}y^{b_{1}}-\frac{1}{\mu}k_{2}x^{a_{2}}y^{b_{2}}\right),\quad\dot{y}=(\widetilde{b}_{1}-b_{1})\left(k_{1}x^{a_{1}}y^{b_{1}}-\frac{1}{\mu}k_{2}x^{a_{2}}y^{b_{2}}\right).

(3.4)

If $(a_{1},b_{1})=(a_{2},b_{2})$ , then

\displaystyle\dot{x}=(\widetilde{a}_{1}-a_{1})\left(k_{1}-\frac{1}{\mu}k_{2}\right)x^{a_{1}}y^{b_{1}},\quad\dot{y}=(\widetilde{b}_{1}-b_{1})\left(k_{1}-\frac{1}{\mu}k_{2}\right)x^{a_{1}}y^{b_{1}}.

If $k_{1}=k_{2}/\mu$ , then every positive point is a steady state and so the system is not static ACR. If $k_{1}\neq k_{2}/\mu$ , then at least one of $\dot{x}$ or $\dot{y}$ is either positive on all of $\mathbb{R}^{2}_{>0}$ or negative on all of $\mathbb{R}^{2}_{>0}$ . But then there is no positive steady state. So $\mathcal{G}$ does not have the capacity for static ACR.

From (3.4), steady states must satisfy the equation

\displaystyle x^{a_{1}-a_{2}}y^{b_{1}-b_{2}}=\frac{k_{2}}{\mu k_{1}}=:k.

(3.5)

Now suppose that $0\not\in\{a_{2}-a_{1},b_{2}-b_{1}\}$ . From (3.5), we see that

(x_{\beta},y_{\beta})=\left((k\beta)^{1/(a_{1}-a_{2})},\beta^{-1/(b_{1}-b_{2})}\right)

is a steady state for every $\beta\in\mathbb{R}_{>0}$ . In particular, two distinct choices of $\beta$ result in distinct $x$ and $y$ components in the two steady states. This implies that neither variable is static ACR for any $k$ . So $\mathcal{G}$ does not have the capacity for static ACR.

Finally, we show that ( ${\it 4}\implies{\it 3}$ ). Assume without loss of generality that $a_{2}-a_{1}\neq 0$ and $b_{2}-b_{1}=0$ . From (3.5), we have that $x=k^{*}\coloneqq k^{\frac{1}{a_{1}-a_{2}}}$ at steady state, which shows that the system is static ACR and $x$ is the static ACR variable. Since this is true for every choice of mass action rate constants, $\mathcal{G}$ is static ACR and $X$ is a static ACR species. Every point on the hyperplane $\mathcal{H}[x,k^{*}]$ is a steady state which shows that $X$ is strong static ACR.

It is clear that the roles of species $X$ and $Y$ are reversed if we assume that $a_{2}=a_{1}$ and $b_{2}\neq b_{1}$ , which proves the claims about the species $Y$ . ∎

Conditions for static ACR in reaction networks with 2 reactions and 2 or fewer species have also been studied in [18].

3.4 Network motifs and their embeddings

Similar to static ACR, we will show that dynamic ACR is a network property. If a network has the capacity for dynamic ACR, then it is dynamic ACR. Moreover, whether a network has the capacity for dynamic ACR depends only on its topology and not on the specific embedding in the Euclidean plane. We refer to such a class of networks as a motif. A network motif in two dimensions is determined by: (i) slope of the reactant polytope, (ii) the quadrant or axis each reaction points along, and (iii) the relative slopes of the two reactions. We demonstrate the relation between a network motif and its multiple embeddings via an example in Figure 12.

(a) Motif

(b) Embedding 1

Figure 12: A network motif (on top) and two of its embeddings (bottom row).

3.5 Dynamic ACR in static ACR networks with 2 reactions & 2 or fewer species

Theorem 3.6 (Dynamic ACR in static ACR networks with 2 reactions & 2 or fewer species).

Let $\mathcal{G}$ be as in (3.1)-(3.2) and suppose that $\dim(\SS)=1$ . The following statements A1-A6 are equivalent:

A1.

$\mathcal{G}$ has the capacity for weak dynamic ACR.
A2.

$\mathcal{G}$ is weak dynamic ACR.
A3.

$\mathcal{G}$ is full basin weak dynamic ACR.
A4.

$\mathcal{G}$ has the capacity for dynamic ACR.
A5.

$\mathcal{G}$ is (full space) dynamic ACR.
A6.

$\mathcal{G}$ is static ACR and $(\widetilde{a}_{1}-a_{1})(a_{2}-a_{1})+(\widetilde{b}_{1}-b_{1})(b_{2}-b_{1})>0$ .

Now, suppose that $\mathcal{G}$ has 2 species and is dynamic ACR. Then either $X$ or $Y$ , but not both, is a dynamic ACR species. Furthermore, the following statements B1-B3 are equivalent.

B1.

$X$ is a dynamic ACR species.
B2.

$X$ is a static ACR species.
B3.

$a_{2}\neq a_{1}$ .

•

When $X$ is a dynamic ACR species,

\mbox{the dynamic ACR value of the variable }x=\mbox{ the static ACR value of }x=\left(\frac{k_{2}}{\mu k_{1}}\right)^{\frac{1}{a_{1}-a_{2}}}.

Analogous statements to B1-B3 and the statement about ACR value hold when $X$ is replaced with $Y$ , and $a_{i}$ is replaced with $b_{i}$ for $i\in\{1,2\}$ .

Proof.

( ${\it A3}\implies{\it A2}\implies{\it A1}$ ) and ( ${\it A5}\implies{\it A4}\implies{\it A1}$ ) hold by definition.

We now show that ( ${\it A1}\implies{\it A6}\implies{\it A5,A3}$ ). Suppose that $\mathcal{G}$ has the capacity for weak dynamic ACR. From properties of one dimensional dynamical systems, it is clear that $\mathcal{G}$ has the capacity for static ACR. From Theorem 3.5, $\mathcal{G}$ is static ACR and one of the two (but not both) is a static ACR species. Assume that $X$ (and not $Y$ ) is the static ACR species. Then $b_{2}=b_{1}$ , $a_{2}\neq a_{1}$ , and $x$ is static ACR with value $\displaystyle x^{*}\coloneqq\left(k_{2}/(\mu k_{1})\right)^{\frac{1}{a_{1}-a_{2}}}$ .

\displaystyle\dot{x}=(\widetilde{a}_{1}-a_{1})y^{b_{1}}x^{a_{1}}\left(k_{1}-\frac{1}{\mu}k_{2}x^{a_{2}-a_{1}}\right).

Clearly, the steady state $x^{*}$ is stable if and only if $(\widetilde{a}_{1}-a_{1})(a_{2}-a_{1})>0$ . If we assume instead that $Y$ (and not $X$ ) is the static ACR species, we get the stability condition $(\widetilde{b}_{1}-b_{1})(b_{2}-b_{1})>0$ . The desired inequality in A6 is obtained by combining the two stability conditions, since it is always the case that one term in A6 is positive and the other term is zero which shows that ( ${\it A1}\implies{\it A6}$ ). Since a unique (within a compatibility class) steady state that is stable must be globally stable for a one dimensional system (i.e. attracts all compatible positive points), we also have that ( ${\it A6}\implies{\it A5}$ ). Moreover, the initial values that are not compatible with the hyperplane of steady states $\{x=x^{*}\}$ also result in trajectories that move towards $\{x=x^{*}\}$ but converge at a boundary steady state. This gives full basin weak dynamic ACR, so we have also proved that ( ${\it A6}\implies{\it A3}$ ).

The last part also shows that when $\mathcal{G}$ is dynamic ACR, (B3 $\implies$ B1 $\implies$ B2) and (B2 $\implies$ B3) is from Theorem 3.5. The statement about the dynamic ACR value of $x$ also follows from the last part. Clearly, the assumption that $a_{2}=a_{1}$ and $b_{2}\neq b_{1}$ will switch the roles of $X$ and $Y$ and give analogous statement for species $Y$ . ∎

Figure 13: Network motifs with 2 reactions and 1 dimensional stoichiometric subspace that show static ACR. The three network types in the magenta ellipse are dynamic ACR, the two network types in the smaller cyan ellipse are wide basin dynamic ACR, while the one network type in the smallest teal ellipse is full basin dynamic ACR. In each of the networks, the static ACR species is

X

, which is on the axis parallel to the reaction polytope.

Theorem 3.7 (Wide basin dynamic ACR in 2 reaction, 2 or fewer species networks with $\dim(\SS)=1$ ).

Let $\mathcal{G}$ be as in (3.1)-(3.2). Suppose that $\dim(\SS)=1$ and $\mathcal{G}$ is dynamic ACR.

1.

The dynamic ACR species is wide basin if and only if $(\widetilde{a}_{1}-a_{1})(\widetilde{b}_{1}-b_{1})\leq 0$ .
2.

The dynamic ACR species is full basin if and only if $(\widetilde{a}_{1}-a_{1})(\widetilde{b}_{1}-b_{1})=0$ .

Proof.

Since $\mathcal{G}$ is dynamic ACR and $\dim(\SS)=1$ , by Theorem 3.6, $\mathcal{G}$ is static ACR. By Theorem 3.5, $v_{1}=-\mu v_{2}$ for some $\mu>0$ , so that the mass action system can be written as in (3.4). This implies that $(\widetilde{b}_{1}-b_{1})\dot{x}-(\widetilde{a}_{1}-a_{1})\dot{y}=0$ , and so two points $(x^{*},y^{*})$ and $(x_{0},y_{0})$ are compatible if and only if

		$\displaystyle(\widetilde{b}_{1}-b_{1})x^{}-(\widetilde{a}_{1}-a_{1})y^{}=(\widetilde{b}_{1}-b_{1})x_{0}-(\widetilde{a}_{1}-a_{1})y_{0}$
		$\displaystyle\iff(\widetilde{b}_{1}-b_{1})(x_{0}-x^{})=(\widetilde{a}_{1}-a_{1})(y_{0}-y^{})$		(3.6)

Since at least one of $\widetilde{a}_{1}-a_{1}$ or $\widetilde{b}_{1}-b_{1}$ must be nonzero, the above implies that $\operatorname{sgn}((x_{0}-x^{*})(y_{0}-y^{*}))=\operatorname{sgn}((\widetilde{b}_{1}-b_{1})(\widetilde{a}_{1}-a_{1}))$ , where $\operatorname{sgn}$ is the sign function that has range $\{+,-,0\}$ and is defined via:

\displaystyle\operatorname{sgn}(z)=\begin{cases}+&\mbox{ if }z>0\\ -&\mbox{ if }z<0\\ 0&\mbox{ if }z=0\end{cases}

We will assume without loss of generality that $X$ is the dynamic ACR species, so that for a particular choice of $K$ , $x$ is the dynamic ACR variable with ACR value $k^{*}$ that depends on $K$ .

Suppose that $\operatorname{sgn}((\widetilde{b}_{1}-b_{1})(\widetilde{a}_{1}-a_{1}))=+$ . Then $\alpha\coloneqq(\widetilde{a}_{1}-a_{1})/(\widetilde{b}_{1}-b_{1})>0$ . From (3.5), $(x_{0},y_{0})$ is compatible with $(x^{*},y^{*})$ if and only if $x_{0}-x^{*}=\alpha(y_{0}-y^{*})$ with $\alpha>0$ . Any initial value $(x_{0},y_{0})$ with $x_{0}>k^{*}+\alpha y_{0}$ is incompatible with $\{x=k^{*}\}$ . To see this, note that

x_{0}-k^{*}>\alpha y_{0}>\alpha(y_{0}-y^{*}).

Clearly $x_{0}$ is not bounded above on the incompatible set $\{x_{0}>k^{*}+\alpha y_{0}\}$ , so in this case $x$ is a narrow basin dynamic ACR variable for any choice of $K$ . In other words, $X$ is a narrow basin dynamic ACR species.

Now suppose that $\operatorname{sgn}((\widetilde{b}_{1}-b_{1})(\widetilde{a}_{1}-a_{1}))=0$ . Since $x$ is the dynamic ACR variable, $b_{2}=b_{1}$ . By condition in the previous part, $\widetilde{a}_{1}\neq a_{1}$ and so $\widetilde{b}_{1}=b_{1}$ . In particular $\dot{y}=0$ and every positive initial value is compatible with the set $\{x=k^{*}\}$ . So $x$ is a full basin dynamic ACR variable for every $K$ . In other words, $X$ is a full basin dynamic ACR species.

Finally, suppose that $\operatorname{sgn}((\widetilde{b}_{1}-b_{1})(\widetilde{a}_{1}-a_{1}))=-$ . In this case, $x_{0}-x^{*}=-\beta(y_{0}-y^{*})$ with $\beta=-(\widetilde{a}_{1}-a_{1})/(\widetilde{b}_{1}-b_{1})>0$ . Any initial value $(x_{0},y_{0})$ with $x_{0}+\beta y_{0}>k^{*}$ is compatible with $\{x=k^{*}\}$ . To see this, let

y^{*}\coloneqq(x_{0}+\beta y_{0}-k^{*})/\beta>0.

Then $(x_{0},y_{0})$ is compatible with $(k^{*},y^{*})$ . So $x$ is a wide basin dynamic ACR variable for every $K$ . In other words, $X$ is a wide basin dynamic ACR species. The only remaining thing to show is that $X$ is not full basin. Let $(x_{0},y_{0})\coloneqq\left(k^{*}/4,k^{*}/(4\beta)\right)\in\mathbb{R}^{2}_{>0}$ . Then

x_{0}+\beta y_{0}=k^{*}/2<k^{*}<k^{*}+\beta y^{*},

and so $(x_{0},y_{0})$ is not compatible with $\{x=k^{*}\}$ .

This completes the proof since we have covered the entire range of $\operatorname{sgn}((\widetilde{b}_{1}-b_{1})(\widetilde{a}_{1}-a_{1}))$ . ∎

The results of Theorems 3.5–3.7 can be translated into a pictorial representation using the idea of motifs (see Figure 13). There are 8 distinct motifs, each has two reactions depicted with red arrows. The reactants are connected by a green line segment, the reactant polytope of the reaction network. See [16] where this representation for networks with two reactions was used to identify the multistationarity property. An upward pointing arrow indicates a reaction of the type $mY\to nY$ for some $n>m$ . An arrow pointing towards the north-west indicates a reaction of the type $a_{1}X+b_{1}Y\to a_{2}X+b_{2}Y$ with $a_{2}<a_{1}$ and $b_{2}>b_{1}$ , and so on.

4 Dynamic ACR in networks with an invariant hyperplane

We now consider networks with 2 reactions and 2 or fewer species and allow the stoichiometric space to have dimension either 1 or 2. With dimension 2 or higher, there is a possibility of oscillatory solutions. However, it is easy to see that at least three reactions are required for such solutions (see also [21]). It is worth mentioning in passing that for general reaction networks, oscillations are compatible with both static and dynamic ACR. An example of a mass action system with three reactions that has oscillations and static ACR in both species is the classic Lotka-Volterra system ( $A\to 2A,A+B\to 2B,B\to 0$ ), see [13] for a discussion on this. The same paper [13] also features an example with oscillations and dynamic ACR, but it requires 6 reactions and a three dimensional stoichiometric subspace. When a network with 2 reactions and 2 species has static ACR, every point on the hyperplane $\mathcal{H}[x,x^{*}]$ is a steady state. A natural generalization of this property is to consider the family of networks for which the hyperplane $\mathcal{H}[x,x^{*}]$ is invariant. For $\dim(\SS)=2$ networks, static ACR is ruled out by results from the previous section but dynamic ACR is still possible.

Theorem 4.1.

Let $(\mathcal{G},K)$ be as in (3.1)-(3.2).

1.

There is a unique $x^{*}>0$ such that $\mathcal{H}[x,x^{*}]$ is invariant if and only if $b_{2}=b_{1}$ , $a_{2}\neq a_{1}$ and $(\widetilde{a}_{1}-a_{1})(\widetilde{a}_{2}-a_{2})<0$ .
2.

There is a unique $x^{*}>0$ such that $\mathcal{H}[x,x^{*}]$ is globally weakly attracting if and only if $\mathcal{H}[x,x^{*}]$ is globally weakly stable if and only if $b_{2}=b_{1}$ , $a_{2}\neq a_{1}$ , $(\widetilde{a}_{1}-a_{1})(\widetilde{a}_{2}-a_{2})<0$ and $(a_{2}-a_{1})(\widetilde{a}_{1}-a_{1})>0$ .

Proof.

The condition that $0=\dot{x}|_{x=x^{*}}$ is equivalent to:

\displaystyle 0=k_{1}(\widetilde{a}_{1}-a_{1})+k_{2}(\widetilde{a}_{2}-a_{2})(x^{*})^{a_{2}-a_{1}}y^{b_{2}-b_{1}}.

(4.1)

For any $x=x^{*}$ , the equation is an identity if $\widetilde{a}_{1}=a_{1}$ and $\widetilde{a}_{2}=a_{2}$ . So assume that this is not the case. It is clear that there is a unique positive $x^{*}$ for which (4.1) holds if and only if $a_{2}\neq a_{1}$ , $b_{2}=b_{1}$ and $(\widetilde{a}_{1}-a_{1})(\widetilde{a}_{2}-a_{2})<0$ .

For the second part, it suffices to assume that $a_{2}\neq a_{1}$ , $b_{2}=b_{1}$ and $(\widetilde{a}_{1}-a_{1})(\widetilde{a}_{2}-a_{2})<0$ because otherwise $\mathcal{H}[x,x^{*}]$ is not invariant and so is not globally weakly attracting. Further we may assume that $a_{2}>a_{1}$ , possibly after relabeling of reactions. With these assumptions, the original dynamical system $\mathcal{D}$ in (3.2) has the same trajectories as the following dynamical system $\mathcal{D}^{\prime}$ :

	$\displaystyle\dot{x}$	$\displaystyle=k_{1}(\widetilde{a}_{1}-a_{1})+k_{2}(\widetilde{a}_{2}-a_{2})x^{a_{2}-a_{1}}$
	$\displaystyle\dot{y}$	$\displaystyle=k_{1}(\widetilde{b}_{1}-b_{1})+k_{2}(\widetilde{b}_{2}-b_{2})x^{a_{2}-a_{1}}.$		(4.2)

Note that the $x$ equation is autonomous and does not have a $y$ -dependence. So a solution to the $x$ equation exists for all time $t\in[0,\infty)$ . If $\widetilde{a}_{1}<a_{1}$ and $\widetilde{a}_{2}>a_{2}$ then the $x$ -component of the trajectories moves away from $x^{*}$ monotonically and so $x^{*}$ is not weakly attracting. If $\widetilde{a}_{1}>a_{1}$ and $\widetilde{a}_{2}<a_{2}$ then the $x$ -component of the trajectories moves towards $x^{*}$ . This proves the second part. ∎

Theorem 4.2.

Let $(\mathcal{G},K)$ be as in (3.1)-(3.2). Suppose that $\mathcal{H}[x,x^{*}]$ is globally weakly attracting for some $x^{*}>0$ . Denote by $\sigma_{i}$ the slope of the reaction vector $v_{i}$ , $i\in\{1,2\}$ .

1.

$x(t)\not\to\mathcal{H}[x,x^{*}]$ for any $(x(0),y(0))\not\in\mathcal{H}[x,x^{*}]$ if $(a_{2}-a_{1})(\sigma_{2}-\sigma_{1})<0$ .
2.

$x(t)\to\mathcal{H}[x,x^{*}]$ for every $(x(0),y(0))\in\mathcal{H}[x,x^{*}]+{\rm span}\{v_{2}\}$ if $(a_{2}-a_{1})(\sigma_{2}-\sigma_{1})\geq 0$ .
3.

$x(t)\to\mathcal{H}[x,x^{*}]$ for every initial value in some cylindrical neighborhood of $\mathcal{H}[x,x^{*}]$ if $(a_{2}-a_{1})(\sigma_{2}-\sigma_{1})>0$ .
4.

$\mathcal{H}[x,x^{*}]$ is globally attracting if and only if $\widetilde{b}_{1}\geq b_{1}$ and $\widetilde{b}_{2}\geq b_{2}$ .

Proof.

From the proof of Theorem 4.1, the dynamical system $\mathcal{D}$ is equivalent to $\mathcal{D}^{\prime}$ in (4) with $a_{2}>a_{1}$ , $\widetilde{a}_{1}>a_{1}$ and $\widetilde{a}_{2}<a_{2}$ . Note that for $\mathcal{D}^{\prime}$ , $x\to x^{*}$ and $\dot{x}\to 0$ asymptotically. In particular, $\absolutevalue{\dot{y}}$ is bounded for every initial value. This means that the solution to the dynamical system $\mathcal{D}^{\prime}$ exists for all nonnegative times for every initial value in $\mathbb{R}^{2}_{\geq 0}$ .

Now, the $\dot{y}$ equation may be written as

\displaystyle\dot{y}=k_{1}(\widetilde{a}_{1}-a_{1})(\sigma_{2}-\sigma_{1})+\sigma_{2}\dot{x}

(4.3)

If $\sigma_{2}=\sigma_{1}$ , then every trajectory converges to some steady state, and this positive steady state is on $\mathcal{H}[x,x^{*}]$ when the initial value is compatible with $\mathcal{H}[x,x^{*}]$ . If $\sigma_{2}\neq\sigma_{1}$ , then after some finite time $\dot{y}$ is bounded away from $0$ and $\operatorname{sgn}(\dot{y})=\operatorname{sgn}(\sigma_{2}-\sigma_{1})$ . In particular, for $\sigma_{2}<\sigma_{1}$ , $y(t)$ reaches $0$ in finite time. Therefore, the trajectory fails to converge to $\mathcal{H}[x,x^{*}]$ for any positive initial value not on $\mathcal{H}[x,x^{*}]$ . On the other hand, if $\sigma_{2}>\sigma_{1}$ , then $\dot{y}>\sigma_{2}\dot{x}$ everywhere in $\mathbb{R}^{2}_{>0}$ . When $\dot{y}=\sigma_{2}\dot{x}$ , we have convergence to $\mathcal{H}[x,x^{*}]$ for every $(x(0),y(0))\in\mathcal{H}[x,x^{*}]+{\rm span}\{v_{2}\}$ . So it follows that when $\dot{y}>\sigma_{2}\dot{x}$ , after a finite time, the trajectory enters a cylinder $\mathcal{C}_{\varepsilon}$ with $\varepsilon$ such that $\dot{y}>k_{1}(\widetilde{a}_{1}-a_{1})(\sigma_{2}-\sigma_{1})/2$ within $\mathcal{C}_{\varepsilon}$ . Moreover, every trajectory with an initial value in $\mathcal{C}_{\varepsilon}$ must converge to the hyperplane $\mathcal{H}[x,x^{*}]$ . This completes the proof of points 1-3.

If $\widetilde{b}_{1}\geq b_{1}$ and $\widetilde{b}_{2}\geq b_{2}$ , then it is immediate from (4) that $\dot{y}\geq 0$ everywhere in $\mathbb{R}^{2}_{\geq 0}$ and so $\mathcal{H}[x,x^{*}]$ is globally attracting. To show the converse, suppose first that $\widetilde{b}_{1}<b_{1}$ . There is a $\delta>0$ such that for $\{(x(0),y(0))\in\mathbb{R}^{2}_{>0}:x(0)^{2}+y(0)^{2}<\delta^{2}\}$ , $\dot{y}$ is negative and bounded away from zero. So such a trajectory reaches the $y=0$ boundary and fails to converge to $\mathcal{H}[x,x^{*}]$ . A similar argument works for $\widetilde{b}_{2}<b_{2}$ , where we can choose an initial $x(0)$ sufficiently large and an initial $y(0)$ sufficiently small. Thus in either case $\mathcal{H}[x,x^{*}]$ is not globally attracting. ∎

Theorem 4.3.

Let $\mathcal{G}$ be as in (3.1)-(3.2). Suppose that $X$ is subspace dynamic ACR but not full basin dynamic ACR. Then the following hold:

1.

$X$ is a wide basin dynamic ACR species if and only if $(\widetilde{b}_{i}-b_{i})(a_{i}-a_{j})>0$ for $i\neq j\in\{1,2\}$ .
2.

$X$ is a narrow basin dynamic ACR species if and only if $(\widetilde{b}_{i}-b_{i})(a_{i}-a_{j})<0$ for $i\neq j\in\{1,2\}$ .

Proof.

Since $X$ is subspace dynamic ACR, $(a_{2}-a_{1})(\sigma_{2}-\sigma_{1})\geq 0$ by the previous theorem. We may assume that $a_{2}>a_{1}$ , $\sigma_{2}\geq\sigma_{1}$ , $\widetilde{a}_{1}>a_{1}$ and $\widetilde{a}_{2}<a_{2}$ . The condition in the first statement is equivalent to $\widetilde{b}_{1}<b_{1}$ and $\widetilde{b}_{2}>b_{2}$ , i.e. $\sigma_{2}<0$ . This means that the set of points incompatible with $\mathcal{H}[x,x^{*}]$ is to the left of the hyperplane and therefore the condition is equivalent to wide basin dynamic ACR. A similar argument in the second case shows equivalence with $\sigma_{2}>0$ and therefore with narrow basin dynamic ACR. ∎

5 Summary of ACR properties for networks with 2 reactions and at most 2 species

The following theorem is a summary of the main results on the different ACR properties in networks with 2 reactions and at most 2 species.

Theorem 5.1.

Let $\mathcal{G}$ be a reaction network with 2 species $\{X,Y\}$ and the following two reactions:

\displaystyle a_{1}X+b_{1}Y\xrightarrow{k_{1}}\widetilde{a}_{1}X+\widetilde{b}_{1}Y,\quad a_{2}X+b_{2}Y\xrightarrow{k_{2}}\widetilde{a}_{2}X+\widetilde{b}_{2}Y,

(5.1)

where $a_{i},b_{i},\widetilde{a}_{i},\widetilde{b}_{i}\in\mathbb{R}_{\geq 0}$ and $(\widetilde{a}_{i},\widetilde{b}_{i})\neq(a_{i},b_{i})$ for $i\in\{1,2\}$ . The labels $k_{1}$ and $k_{2}$ are mass action reaction rate constants, when considering a mass action system $(\mathcal{G},K)$ .

1.

A necessary and sufficient condition for the existence of a unique invariant hyperplane $\mathcal{H}$ parallel to a coordinate axis is that the reactant polytope be a line segment parallel to the other coordinate axis. In particular, $(a_{1},b_{1})\neq(a_{2},b_{2})$ and $(a_{2}-a_{1})(b_{2}-b_{1})=0$ .

\suspend
enumerate

Suppose that $\mathcal{G}$ has a unique invariant hyperplane $\mathcal{H}$ parallel to a coordinate axis. Then the following hold:

\resume
enumerate
2.

$\mathcal{G}$ has the capacity for $\mathcal{P}$ -ACR if and only if $\mathcal{G}$ is $\mathcal{P}$ -ACR, where $\mathcal{P}\in\{$ static, strong static, dynamic, weak dynamic $\}$ , i.e. the $\mathcal{P}$ -ACR properties are independent of choice of rate constants.

We use the convention $0\cdot(1/0)=0$ in the following expressions.

(a)

$\mathcal{G}$ is static ACR if and only if $\mathcal{G}$ is strong static ACR if and only if $\dim(\SS)=1$ and reaction vectors are negative scalar multiples of each other: $\widetilde{a}_{1}-a_{1}=-\mu\left(\widetilde{a}_{2}-a_{2}\right)\mbox{ for some }\mu>0.$
(b)

$\mathcal{H}$ is weakly attracting if and only if both reactions point inwards:

$\displaystyle(\widetilde{a}_{i}-a_{i})(a_{j}-a_{i})+(\widetilde{b}_{i}-b_{i})(b_{j}-b_{i})>0,\quad i\neq j\in\{1,2\}.$

(c)

$\mathcal{G}$ is non-null dynamic ACR if and only if $\mathcal{G}$ is subspace dynamic ACR if and only if

\displaystyle\sum_{i\neq j\in\{1,2\}}\left\{(a_{j}-a_{i})\left(\frac{\widetilde{b}_{i}-b_{i}}{\widetilde{a}_{i}-a_{i}}\right)+(b_{j}-b_{i})\left(\frac{\widetilde{a}_{i}-a_{i}}{\widetilde{b}_{i}-b_{i}}\right)\right\}\geq 0.

\suspend

enumerate

Suppose that neither $\dot{x}$ nor $\dot{y}$ is identically zero, or equivalently $(\widetilde{a}_{1},\widetilde{a}_{2})\neq(a_{1},a_{2})$ and $(\widetilde{b}_{1},\widetilde{b}_{2})\neq(b_{1},b_{2})$ .

\resume

enumerate

(d)

$\mathcal{G}$ is cylinder dynamic ACR if and only if

\displaystyle\sum_{i\neq j\in\{1,2\}}\left\{(a_{j}-a_{i})\left(\frac{\widetilde{b}_{i}-b_{i}}{\widetilde{a}_{i}-a_{i}}\right)+(b_{j}-b_{i})\left(\frac{\widetilde{a}_{i}-a_{i}}{\widetilde{b}_{i}-b_{i}}\right)\right\}>0.

(e)

$\mathcal{G}$ is full basin dynamic ACR if and only if

\displaystyle(a_{j}-a_{i})\left(\frac{\widetilde{b}_{i}-b_{i}}{\widetilde{a}_{i}-a_{i}}\right)+(b_{j}-b_{i})\left(\frac{\widetilde{a}_{i}-a_{i}}{\widetilde{b}_{i}-b_{i}}\right)\geq 0,\quad i\neq j\in\{1,2\}.

4.
For all the above networks with $\mathcal{P}$ -ACR:
- •
  
  Either $X$ or $Y$ , but not both, is an $\mathcal{P}$ -ACR species.
- •
  
  $X$ is an $\mathcal{P}$ -ACR species if $a_{2}\neq a_{1}$ . The variable $x$ has the $\mathcal{P}$ -ACR value $\displaystyle\left(-\frac{k_{2}(\widetilde{a}_{2}-a_{2})}{k_{1}(\widetilde{a}_{1}-a_{1})}\right)^{\frac{1}{a_{1}-a_{2}}}$ .
- •
  
  $Y$ is an $\mathcal{P}$ -ACR species if $b_{2}\neq b_{1}$ . The variable $y$ has the $\mathcal{P}$ -ACR value $\displaystyle\left(-\frac{k_{2}(\widetilde{b}_{2}-b_{2})}{k_{1}(\widetilde{b}_{1}-b_{1})}\right)^{\frac{1}{b_{1}-b_{2}}}$ .

A pictorial representation of these results is shown in Figure 14.

6 Discussion and future work

In this paper, we have established that for small networks, the Euclidean embedding (or geometric structure) of a reaction network can yield deep insights into the dynamics of the mass action ODE system. In particular, when there are only two reactions and at most two species, the reactant polytope (the line segment joining the two reactant complexes) is required to be horizontal or vertical for any type of ACR property. Some of the two reaction networks that appear in this paper do show up in applications. For instance, the archetypal wide basin ACR network in 2(a) can be thought of as a simple model of infectious disease dynamics (SIS model). The reaction $A+B\to 2B$ represents an infective individual $B$ infecting a susceptible individual $A$ , while the reaction $B\to A$ represents recovery from infection. One may also interpret this model as a protein with two alternate conformations $A$ and $B$ with the spontaneous transition $B\to A$ , while the transition $A\to B$ is catalyzed (or promoted) by $B$ . However, most biochemically realistic networks are significantly more complex with several reactions and species. Even though ACR can be found in higher dimensional systems, we do not expect that there will be such simple characterization of their ACR properties. However, our work suggests that it may be fruitful to study the link between the geometry of reactant polytopes and dynamics of ACR systems further. Moreover, we do expect that small motifs that are embedded within large and complicated networks may affect the overall dynamics. Our future goal is to understand such effects.

The small motifs studied in this paper also serve as test cases for various dynamical behaviors. For instance, a surprising possibility revealed from the study of small motifs was that of weak dynamic ACR, where every trajectory monotonically approaches a hyperplane while simultaneously failing to converge (see Figure 3(b)).

In future work, we study dynamic ACR in more complex, biochemically realistic systems such as bacterial two-component signaling systems, a class that encompasses several thousands of systems [1, 6, 2, 7, 24, 25]. Moreover, we study consequences of dynamic ACR. We plan to show that dynamic ACR, i.e. the property of robustness against variations in initial conditions, surprisingly leads to other much more robust dynamical properties with stronger implications for biochemical systems with dynamic ACR.

Figure 14: (Motifs of Weak Dynamic ACR) There are 17 motifs of weak dynamic ACR with two reactions and two or fewer species. A necessary and sufficient condition for weak dynamic ACR (as well as weak full basin dynamic ACR) is that the reactant polytope be parallel to the axis of the ACR variable ( green line segment) and both reactions point inwards. Of these motifs, 16 are placed on the circumference of a circle with coordinates

\theta=n\pi/8,n\in\{0,1,\ldots 15\}

, while 1 motif is placed at the center of the circle. The two arrows make the same angle with the reactant polytope for the motifs at

\theta=n\pi/2,n\in\{0,1,2,3\}

(the four cardinal directions – north, south, east and west). For

\theta\in n\pi/8,n\in\{0,1,2,3,4\}

(north-east quadrant of the picture), the left arrow is fixed in the north-east quadrant while the right arrow rotates southwards moving southwards along the picture. Similarly in the north-west quadrant, the right arrow is fixed in the north-west quadrant; in the south-west quadrant the left arrow is fixed in the south-east quadrant; and in the south-east quadrant the right arrow is fixed in the south-west quadrant. The figure of motifs is invariant under reflection and under rotation around the central vertical axis. The figure is also invariant under a combination of reflection around a central horizontal axis and rotation of each motif around the axis of that motif. The central horizontal band – with the motifs at

\theta=0,\pi

on the circumference and the motif at the center of the circle – have

\dim(\SS)=1

while the rest have

\dim(\SS)=2

. The motif at the center of the circle can have an embedding with either 1 or 2 species (the second species remaining dynamically unchanged), an embedding of every motif on the circumference requires 2 species. Each motif is labeled with its strongest local ACR property (in magenta) as well as its strongest non-local ACR property (in cyan). Null ACR is labeled in (in olive). Moving northwards along the circumference of the circle, the motifs have stronger local and non-local ACR properties.

Acknowledgments

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

References

[1] Uri Alon. An introduction to systems biology: design principles of biological circuits. CRC press, 2019.
[2] Uri Alon, Michael G Surette, Naama Barkai, and Stanislas Leibler. Robustness in bacterial chemotaxis. Nature, 397(6715):168–171, 1999.
[3] David F Anderson, Daniele Cappelletti, and Thomas G Kurtz. Finite time distributions of stochastically modeled chemical systems with absolute concentration robustness. SIAM Journal on Applied Dynamical Systems, 16(3):1309–1339, 2017.
[4] David F Anderson, Germán A Enciso, and Matthew D Johnston. Stochastic analysis of biochemical reaction networks with absolute concentration robustness. Journal of The Royal Society Interface, 11(93):20130943, 2014.
[5] Murad Banaji. Inheritance of oscillation in chemical reaction networks. Applied Mathematics and Computation, 325:191–209, 2018.
[6] Naama Barkai and Stan Leibler. Robustness in simple biochemical networks. Nature, 387(6636):913–917, 1997.
[7] 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.
[8] 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.
[9] Gheorghe Craciun, Badal Joshi, Casian Pantea, and Ike Tan. Multistationarity in cyclic sequestration-transmutation networks. arXiv preprint arXiv:2110.13975, 2021.
[10] Joseph P Dexter, Tathagata Dasgupta, and Jeremy Gunawardena. Invariants reveal multiple forms of robustness in bifunctional enzyme systems. Integrative Biology, 7(8):883–894, 2015.
[11] Joseph P Dexter and Jeremy Gunawardena. Dimerization and bifunctionality confer robustness to the isocitrate dehydrogenase regulatory system in escherichia coli. Journal of Biological Chemistry, 288(8):5770–5778, 2013.
[12] Badal Joshi. Complete characterization by multistationarity of fully open networks with one non-flow reaction. Applied Mathematics and Computation, 219:6931–6945, 2013.
[13] Badal Joshi and Gheorghe Craciun. Foundations of Static and Dynamic Absolute Concentration Robustness. arXiv preprint arXiv:2104.14070, 2021.
[14] Badal Joshi and Anne Shiu. Atoms of multistationarity in chemical reaction networks. Journal of Mathematical Chemistry, 51(1):153–178, 2013.
[15] 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.
[16] Badal Joshi and Anne Shiu. Which small reaction networks are multistationary? SIAM Journal on Applied Dynamical Systems, 16(2):802–833, 2017.
[17] Robert L Karp, Mercedes Pérez Millán, Tathagata Dasgupta, Alicia Dickenstein, and Jeremy Gunawardena. Complex-linear invariants of biochemical networks. Journal of theoretical biology, 311:130–138, 2012.
[18] Nicolette Meshkat, Anne Shiu, and Angelica Torres. Absolute concentration robustness in networks with low-dimensional stoichiometric subspace. Vietnam Journal of Mathematics, pages 1–29, 2021.
[19] B Pascual-Escudero and E Feliu. Local and global robustness in systems of polynomial equations. arXiv preprint arXiv:2005.08796, 2020.
[20] Friedrich Schlögl. Chemical reaction models for non-equilibrium phase transitions. Zeitschrift für Physik, 253(2):147–161, 1972.
[21] J Schnakenberg. Simple chemical reaction systems with limit cycle behaviour. Journal of theoretical biology, 81(3):389–400, 1979.
[22] Guy Shinar and Martin Feinberg. Structural sources of robustness in biochemical reaction networks. Science, 327(5971):1389–1391, 2010.
[23] Guy Shinar and Martin Feinberg. Design principles for robust biochemical reaction networks: what works, what cannot work, and what might almost work. Mathematical biosciences, 231(1):39–48, 2011.
[24] 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.
[25] Guy Shinar, Joshua D Rabinowitz, and Uri Alon. Robustness in glyoxylate bypass regulation. PLoS Comput Biol, 5(3):e1000297, 2009.
[26] Polly Y Yu and Gheorghe Craciun. Mathematical analysis of chemical reaction systems. Israel Journal of Chemistry, 58(6-7):733–741, 2018.

Reaction Network Motifs for Static and Dynamic Absolute Concentration Robustness

Abstract

1 Introduction

2 Forms of ACR

Background and terminology

Definition 2.1.

Definition 2.2.

2.1 Static forms of ACR

Definition 2.3.

Remark 2.4.

Definition 2.5.

2.2 Dynamic forms of ACR

Definition 2.6.

Definition 2.7.

Remark 2.8.

Definition 2.9.

Remark 2.10.

Example 2.11.

2.3 Relations between different forms of ACR

Theorem 2.12.

Proof.

2.4 Some basins of interest for all forms of ACR

Definition 2.13.

Theorem 2.14.

Proof.

2.5 Some examples of systems with local dynamic ACR

2.6 An example of a system with subspace (but not full space) dynamic ACR

3 Classification of Minimal Static and Dynamic ACR Networks

Definition 3.1.

3.1 Static and dynamic ACR for reaction networks with only 1 reaction or only 1 species.

Theorem 3.2 (Static and Dynamic ACR in one-reaction networks).

Proof.

Theorem 3.3 (Static ACR in one-species networks).

Proof.

Theorem 3.4 (Dynamic ACR in one-species networks).

Proof.

3.2 Reaction networks with 2 reactions and 2 or fewer species: Notation

3.3 Static ACR for reaction networks with 2 reactions and 2 or fewer species.

Theorem 3.5 (Static ACR in networks with 2 reactions and 2 or fewer species).

Proof.

3.4 Network motifs and their embeddings

3.5 Dynamic ACR in static ACR networks with 2 reactions & 2 or fewer species

Theorem 3.6 (Dynamic ACR in static ACR networks with 2 reactions & 2 or fewer species).

Proof.

Theorem 3.7 (Wide basin dynamic ACR in 2 reaction, 2 or fewer species networks with dim(SS)=1\dim(\SS)=1).

Proof.

4 Dynamic ACR in networks with an invariant hyperplane

Theorem 4.1.

Proof.

Theorem 4.2.

Proof.

Theorem 4.3.

Proof.

5 Summary of ACR properties for networks with 2 reactions and at most 2 species

Theorem 5.1.

6 Discussion and future work

Acknowledgments

References

Theorem 3.7 (Wide basin dynamic ACR in 2 reaction, 2 or fewer species networks with $\dim(\SS)=1$ ).