Graphical characterizations of robust stability in biological interaction networks

Abstract

Previous studies have inferred robust stability of reaction networks by utilizing linear programs or iterative algorithms. Such algorithms become tedious or computationally infeasible for large networks. In addition, they operate like black boxes without offering intuition for the structures that are necessary to maintain stability. In this work, we provide several graphical criteria for constructing robust stability certificates, checking robust non-degeneracy, verifying persistence, and establishing global stability. By characterizing a set of stability-preserving graph modifications that includes the enzymatic modification motif, we show that the stability of arbitrarily large nonlinear networks can be examined by simple visual inspection. We show applications of this technique to ubiquitous motifs in systems biology such as post-translational modification (PTM) cycles, the ribosome flow model (RFM), T-cell kinetic proofreading, and others. The results of this paper are dedicated in honor of Eduardo D. Sontag’s seventieth birthday and his pioneering work in nonlinear dynamical systems and mathematical systems biology.

1 Introduction

Biomolecular interaction networks (BINs) function under severe forms of external and internal uncertainty. Nevertheless, they operate robustly and consistently to maintain homeostasis, which is understood as the maintenance of a desired steady state against environmental factors, external signals, and in vivo fluctuations in the concentrations of biochemical species. In fact, robustness has been proposed as a key defining property of biological networks [1, 2]. However, mathematical analysis of such networks has been lagging as the dynamical system descriptions of such networks suffer from nonlinearity and uncertainty. Generic nonlinear dynamical systems are already difficult to analyze due to the scarcity of general and powerful analysis tools. Furthermore, they can manifest complex forms of unstable behavior that are not exhibited by linear systems. For instance, small fluctuations in concentrations, or tiny changes in kinetic parameters, can have radical effects causing the observable phenotype to be driven to a different region of the state space, and/or to lose stability altogether and transform into a sustained oscillation or chaotic behavior. This may make the biological network lose its function and cause key species to reach undesirable or even unsafe levels. In fact, disease can be often characterized mathematically as the loss of stability of a certain phenotype [3, 4]. A second complicating factor is the fact that the exact form of kinetics (determining the speed of interactions) is difficult to measure and is subject to environmental changes.

Fig. 1

Computational RLF construction is tedious for large networks. The figure depicts a Petri-net representation of the PTM star: a substrate that is a target of an arbitrary finite number of distinct competing PTM cycles (e.g., phosphorylation, methylation, ubiquitination, etc). The subnetwork inside the dotted triangle depicts a single PTM cycle. A rectangle denotes a reaction, while a circle denotes a species

Therefore, verifying the stability of a given nonlinear BIN without reference to its kinetics has been a challenging long-standing goal in systems biology research [5]. Nevertheless, partial success has been achieved in this endeavor. Examples include the theory of complex balance [6,7,8], and the theory of monotone BINs [9]. More recently, stability certificates have been constructed via robust Lyapunov functions (RLFs) in reaction [10,11,12], and concentration coordinates [12,13,14,15,16]. Except for a small subclass of BINs (see §III.C), such methods mainly utilize computational algorithms to construct RLFs via either iterative algorithms or linear programs. However, such algorithms act as “black boxes” and are not interpretable in terms of the structural properties of the network’s graph. This has several drawbacks. First, computational algorithms become tedious for larger networks as the number of species and reactions grow. Consider the PTM star depicted in Fig.  whose size grows considerably for large n. Second, “stability-preserving” graph modifications are not well characterized. A simple modification of the BIN graph mandates a rerun of the computational algorithm from scratch. For instance, is the stability of the PTM star preserved if we added inflow/outflow reactions for the substrate ($\emptyset \rightleftharpoons \text{ Substrate }$)?. Third, fundamental “motifs” have been described as the building blocks of BINs [17]. However, a corresponding “modular” theory for RLF construction that utilizes the stability properties of its subnetworks is lacking. For example, the difference between the PTM star (Fig. 1) with n and $n+1$ products is in the addition of an extra PTM cycle. How does the addition of the extra motif affect stability?

The above questions are hard to answer using computational algorithms. In this work, we identify a set of stability-preserving graph modifications. In particular, we show that the stability of many large networks in systems biology can be understood modularly. For the specific network in Fig. 1, we will show that it can be “reduced” to a simple linear network (See Fig.   and Sect. 7.1.2). Hence, it admits a stability certificate for every $n \ge 1$, a result which is not readily achievable using previous results [6,7,8,9, 12, 14]. We will show that the addition of an inflow/outflow reaction to the substrate preserves stability, and that the PTM cycle is a fundamental “stable” motif in a precise manner to be defined.

Our unified framework can be applied to many networks in the literature whose stability was studied individually via various techniques, this includes the T-cell kinetic proofreading network [8], the PTM cycle [18], the all-encompassing processive PTM cycle [19], the ribosome flow model and its variations [20,21,22], and others.

It is worth noting that many of the properties of BINs have already been characterized graphically. This includes complex balance [6, 7], injectivity [23, 24], monotonicity [9, 25], and persistence [26]. Additional studies have tackled graph modifications that preserve various other properties of BINs [27, 28]. Therefore, we complement this literature by characterizing robust stability in graphical terms for classes of BINs for the first time.

The paper proceeds as follows. Section 2 reviews notation and definitions. Section 3 reviews relevant results on linear (monomolecular) networks. In Sect. 3, we list the graph modifications under consideration, and show the existence of RLFs for classes of modified networks. Global stability and robust non-degeneracy are discussed in Sect. 4. Applications are studied in Sect. 5. Proofs are included in the appendix.

2 Background and notation

2.1 Biological interaction networks

Any collection of chemical reactions can be written mathematically using the formalism of biological interaction networks (BINs). Hence, we review the standard definitions and notation [7, 8, 12, 29, 30].

A BIN (also known as a chemical reaction network (CRN)) is a pair ${{\mathscr {N}}}=({{\mathscr {S}}},{{\mathscr {R}}})$, where ${\mathscr {S}}=\{X_1,..,X_n\}$ is the set of species, and ${{\mathscr {R}}}=\{{\textbf{R}}_1,...,{\textbf{R}}_\nu \}$ is the set of reactions. A species is the entity that partakes in or is formed in a chemical interaction. Within the realm of biomolecular networks, a species can be a substrate, a complex, an enzyme, an mRNA molecular, a gene promoter state, etc. A reaction is the transformation of reacting species into product species. Examples include complex formation, binding, unbinding, decay, production, complex formation, etc.

The mathematical structure of BINs can be described by two mathematical substructures: the stoichiometry and the kinetics.

2.1.1 The stoichiometry

The relative gain or loss of molecules of species $X_i$ between the sides of each reaction is the stoichiometry of $X_i$. This is represented by writing a reaction as:

$$\begin{aligned} {\textbf{R}}_j: \quad \sum _{i=1}^n \alpha _{ij} X_i \longrightarrow \sum _{i=1}^n \beta _{ij} X_i, \ j=1,..,\nu , \end{aligned}$$

(1)

where $\alpha _{ij}, \beta _{ij} \ge 0$ are integers known as the stoichiometry coefficients. If a transformation can happen also in the reverse direction, then ${\textbf{R}}_j$ is said to be reversible and its reverse is denoted by ${\textbf{R}}_{-j}$. A reaction can have no reactants or no products (though not simultaneously). The empty side is denoted by $\emptyset$.

If a reaction has a species both as a reactant and as a product (for example, $X+Y \rightarrow X$) then it is called catalytic.

The stoichiometry matrix $\Gamma$ of a given network is an $n \times \nu$ matrix whose (i, j)th entry describes the net gain/loss of the ith species at the jth reaction. Hence, it can be written element wise as:

$[\Gamma ] _{ij} = \beta _{ij}- \alpha _{ij}.$

2.1.2 Kinetics

The set of relationships that determine the speed of transformation of reactant species into product species are known as kinetics. In order to describe such relations, the species need to be quantified. A species $X_i$ is quantified by assigning it a nonnegative real number known as the concentration $x_i {\in {\mathbb {R}}_{\ge 0}^n}$, where ${{\mathbb {R}}}_{\ge 0}^n$ denotes the nonnegative orthant in the n-dimensional Euclidean space. A reaction ${\textbf{R}}_j$ is assigned a single-valued mapping $R_j:{{{\mathbb {R}}}_{\ge 0}^n\rightarrow {{\mathbb {R}}}_{\ge 0}}$ known as the reaction rate. The reaction rate vector is written as $R(x)=[R_1(x),...,R_\nu (x)]^T$.

The most common form of kinetics is known as Mass-Action and it can be written as: $R_j(x)= k_j \prod _{i=1}^n x_i^{\alpha _{ij}}$, where $k_j>0, j=1,..,\nu$ are the kinetic constants. However, this form “is not based on fundamental laws” and is merely “good phenomenology” justified by imagining the reactants as colliding molecules [31]. In biological systems, in particular, other forms of kinetics usually arise when modeling networks involving multiple timescales. This includes Michaelis–Menten, Hill kinetics, etc. Therefore, we do not assume a specific functional form of kinetics. We only assume that the kinetics are monotone. More precisely, the reaction rates $R_j(x), j=1,..,\nu$ satisfy:

• AK1. each reaction varies smoothly with respects to its reactants, i.e., R(x) is $C^1$;

• AK2. a reaction requires all its reactants to occur, i.e., if $\alpha _{ij}>0$, then $x_i=0$ implies $R_j(x)=0$;

• AK3. if a reactant increases, then the reaction rate increases, i.e., ${\partial R_j}/{\partial x_i}(x) \ge 0$ if $\alpha _{ij}>0$ and ${\partial R_j}/{\partial x_i}(x)\equiv 0$ if $\alpha _{ij}=0$. Furthermore, the aforementioned inequality is strict whenever the reactants are strictly positive.

For a given network ${{\mathscr {N}}}$, the set of a reaction rates satisfying the assumptions above is called the admissible kinetics. Furthermore, the assumptions AK1–AK3 translate into a sign-pattern constraint on the Jacobian of R. To formalize this, let ${\mathcal {K}}_{{\mathscr {N}}}$ be defined as ${\mathcal {K}}_{{\mathscr {N}}}=\{V\in {\mathbb {R}}^{\nu \times n} \vert [V]_{ji}>0~\text{ whenever }~ X_i~\text{ is } \text{ a } \text{ reactant } \text{ of }~{\textbf{R}}_j, ~\text{ and }~[V]_{ji}=0~\text{ otherwise }\}$. We think of ${\mathcal {K}}_{{\mathscr {N}}}$ as the set of all possible Jacobian matrices $\partial R/\partial x$ evaluated on the positive orthant ${\mathbb {R}}_+^n$.

2.1.3 Dynamics

We view the concentrations as trajectories in time and write them as $x(t)=[x_1(t),...,x_n(t)]^T$. The temporal evolution of the network is given by the following ordinary differential equation (ODE):

$$\begin{aligned} \dot{x} = \Gamma R(x), \; \; x(0)=x_\circ . \end{aligned}$$

(2)

The positive orthant is forward invariant for (2), i.e., if $x_\circ$ is positive, then the trajectory stays positive for all time $t \ge 0$.

In the biomolecular context, there are usually conserved quantities which do not get created or annihilated during the course of the reaction. This can include total amounts of DNA, enzymes, substrates, ribosomes, etc. Mathematically, a stoichiometric conservation law is a nonnegative vector $d \in {\mathbb {R}}_{\ge 0}^n$ satisfying $d^T \Gamma =0$. If d is positive then the network is called conservative.

The existence of a conservation law implies that $d^Tx(t)\equiv d^T x(0)$. Hence, the positive orthant is partitioned into a foliage of subsets known as stoichiometric classes. For a state vector $x_\circ$, the corresponding class is written as ${\mathscr {C}}_{x_\circ }:=(\{x_\circ \}+ \text{ Im }(\Gamma )) \cap {{\mathbb {R}}}_{\ge 0}^n$, and it is forward invariant. Therefore, all Lyapunov functions and claims of stability are relative to a stoichiometric class. For a conservative network, all stoichiometric classes are compact polyhedral sets, and hence all trajectories are bounded. In addition, this guarantees at least one steady state in each stoichiometric class by applying Brouwer’s fixed point theorem to the associated flow of the dynamical system restricted to the stoichiometric class.

A vector v is called a flux if $\Gamma v=0$. In order to simplify the treatment, we will assume the following about the stoichiometry of the network:

• AS1 There exists a positive flux, i.e., $\exists v \in \ker \Gamma$ such that $v\gg 0$.

• AS2 The network has no catalytic reactions.

Assumption AS1 is necessary for the existence of positive steady states for the corresponding dynamical system (2).

2.2 Graphical representation: Petri-nets

BINs can be represented graphically in several ways. We adopt the Petri-net formalism [32] (also known as the species-reaction graph [33]). A Petri-net is a weighted directed bipartite graph. The vertices consist of the set of species ${{\mathscr {S}}}$ (represented by circles) and the set of reactions ${{\mathscr {R}}}$ (represented by rectangles). An edge with a weight w from $X_i\in {{\mathscr {S}}}$ to ${\textbf{R}}_j\in {{\mathscr {R}}}$ means that $X_i$ is a reactant of ${\textbf{R}}_j$ with stoichiometric coefficient w, while the reverse edge means that $X_i$ is a product of ${\textbf{R}}_j$ with a stoichiometric coefficient w. For a more compact representation, if two reactions are the reverse of each other (e.g., ${\textbf{R}}_j,{\textbf{R}}_{-j}$) then they are represented as a single reaction in the Petri-net with reversible edges. In the formalism of Petri-nets [34], the stoichiometric matrix $\Gamma$ is the incidence matrix of the Petri-net.

For example, the PTM star in Fig. 1 corresponds to the following network:

$$\begin{aligned} S+E_i&\rightleftharpoons C_i \longrightarrow P_i + E_i, \end{aligned}$$

(3)

$$\begin{aligned} P_i+F_i&\rightleftharpoons D_i \longrightarrow S + F_i, \end{aligned}$$

(4)

$i=1,..,n$, where S denotes the substrate and $P_i$ denotes the ith product.

2.3 Robust Lyapunov functions

Following our previous work [11,12,13], a locally Lipschitz function $V: {\mathbb {R}}^n \rightarrow {\mathbb {R}}_{\ge 0}$ is a robust Lyapunov function (RLF) for a given network ${{\mathscr {N}}}$ iff:

• 1. it is positive-definite, i.e., $V(x) \ge 0$ for all x, and $V(x)=0$ iff $\Gamma R(x)=0$, and

• 2. it is non-increasing, i.e., $\dot{V}(x) \le 0$ for all x and all R satisfying $\frac{\partial R}{\partial x} \in {\mathcal {K}}_{{\mathscr {N}}}$.

Since V is not assumed to be continuously differentiable, the derivative above is defined in the sense of Dini as $\dot{V}(x):=\limsup _{h \rightarrow 0^+} ( V(x+h\Gamma R(x))-V(x))/h$ [35]. Existence of an RLF guarantees that the steady state set is Lyapunov stable, and that all V’s level sets are trapping [11, 12, 35]. Global stability can be verified by a LaSalle argument or by establishing robust non-degeneracy of the Jacobian [11, 12, 15]. In this paper, we utilize RLFs that can be written as piecewise linear (PWL) functions in terms of the rates. In [12], it has been shown that they can be converted to PWL RLFs in the concentration coordinates and vice versa. Hence, we will subsequently use the term “PWL RLF” to designate an RLF that is piecewise linear either in the rates or in the concentrations.

3 Linear (monomolecular) networks

3.1 Definition and review

Studying general nonlinear BINs is, predictably, a difficult and open problem. In comparison, assuming linearity simplifies the analysis considerably. In order to get a linear ODE with Mass-Action kinetics, all the reactions have to be monomolecular. In other words, there is only a unique reactant with stoichiometry coefficient 1 for each reaction. The resulting ODE can be studied via standard analysis methods for positive linear systems [36, 37], or as a special case of complex-balanced networks [7]. A weaker notion of linearity is a graphical one where the Petri-net is assumed to be linear, [38], which means that each reaction has a unique reactant and a unique product with the stoichiometry coefficients equal to one. Therefore, nonlinear reaction rates are allowed. It has been long observed that the linearity of the Petri-net is sufficient for analysis, i.e., stability analysis can be performed for general monomolecular networks with monotone kinetics [39]. This generalized class of networks is often known as compartmental networks [40]. Hence, we refer to such networks as linear networks since the corresponding Petri-net is linear. Therefore, we use a graphical notion of linearity and not a kinetic one. The definition is stated formally below:

Definition 1

A given BIN ${{\mathscr {N}}}$ is said to be linear if each reaction can be written as either $X_i \mathop {\longrightarrow }\limits ^{R_{ij}} X_j$, $\emptyset \mathop {\longrightarrow }\limits ^{u_i} X_i$, or $X_i \mathop {\longrightarrow }\limits ^{R_{i}} \emptyset$ for some i, j where $R_{ij},R_{i}:{\mathbb {R}}_{\ge 0} \rightarrow {\mathbb {R}}_{\ge 0}, u_i\ge 0$ are the reaction rates.

Applying the assumptions AK1–AK3, we note that $R_{ij}$ can be any single-valued strictly increasing $C^1$ function that vanishes at the origin.

3.2 Existence of Lyapunov functions: sum-of-currents (SoC) RLF

One of the advantages of studying linear networks is that their stability is well characterized. Indeed, it has been long known [39, 40] that linear networks can be studied using a Lyapunov function of the form:

$$\begin{aligned} V(x) = \Vert \dot{x} \Vert _1= \sum _{i=1}^n \left| \sum _{j \ne i} (R_{ji}(x_j) - R_{ij}(x_i)) + u_i - R_i(x_i) \right| , \end{aligned}$$

(5)

where $u_i \ge 0$ is the inflow to species $X_i$. Note that V is PWL in terms of the rates. We state the following theorem that restates the result in [39] using our terminology:

Theorem 1

Let ${{\mathscr {N}}}$ be a linear BIN with any set of admissible reaction rates $\{R_{ij}(x_i),$ $R_i(x_i),u_i\}_{i,j=1}^n$. Let (2) be the associated ODE. Let V be defined as in (5). Then, V is an RLF for ${{\mathscr {N}}}$.

In order to generalize the result above to classes of nonlinear networks, we will provide a new proof of Theorem 1 in the Appendix based on the techniques used in [11, 12, 22]. The same techniques will be generalized to prove Theorem 4. In [12], we have called (5) a sum-of-currents (SoC) RLF, since it is a sum of the absolute values of the currents $dx_i/dt, i=1,..,n$, which is analogous to the electric current $I=dq/dt$, where q is the electric charge.

3.3 Existence of Lyapunov functions: max–min RLF

For a subclass of linear BINs, another Lyapunov function can be used to establish stability, which is the Max–Min RLF [10, 11]. Define the set-valued function: ${\mathcal {R}}(x)=\{R_{ij}(x),R_i(x),u_i \vert i,j=1,..,n, i\ne j \}$. Then, consider the following function:

$$\begin{aligned} V(x)= \max {\mathcal {R}}(x) - \min {\mathcal {R}}(x), \end{aligned}$$

(6)

Note that V is PWL in terms of the rates. The existence of an RLF of the form (6) can be characterized graphically for general BINs [10, 11]. In order to minimize the notational inconvenience, we assume that ${\textbf{1}}$ is a flux for the network ${{\mathscr {N}}}$. Hence, the result can be stated as follows:

Theorem 2

([10, 11]) Let a BIN ${{\mathscr {N}}}$ be given. Assume that it has a unique positive flux equal to ${\textbf{1}}$ and every species $X_i$ is a reactant to a unique reaction. Then, V as defined in (6) is an RLF for ${{\mathscr {N}}}$.

Remark 1

In order to generalize Theorem 2 to accommodate BINs that admit a unique positive flux $v\gg 0$, the reactions in ${\mathcal {R}}(x)$ can be weighed by the corresponding entry in v [11].

4 Stability-preserving graph modifications

4.1 Definitions

Consider a BIN $({\mathscr {S}},{\mathscr {R}})$ that admits an RLF V. Assume that the network is modified to a new network $(\tilde{{\mathscr {S}}}, \tilde{{\mathscr {R}}})$. We are interested in the existence of an RLF for the new network. To be more concrete, we focus on graph modifications listed in Table . As can be noticed, some of these modification can change a linear network into a nonlinear network. First, we formalize the concept of adding an extra product or reactant to a reaction.

Definition 2

Consider a BIN $({{\mathscr {S}}},{{\mathscr {R}}})$. We say that a reaction $\tilde{\textbf{R}}_j$ is an extension of a reaction ${\textbf{R}}_j \in {{\mathscr {R}}}$ if the following holds for each $X_i \in {{\mathscr {S}}}$: if $X_i \in {{\mathscr {S}}}$ is a reactant of ${\textbf{R}}_j \in {{\mathscr {R}}}$, then $X_i$ is a reactant of $\tilde{\textbf{R}}_j \in \tilde{{\mathscr {R}}}$ with the same stoichiometric coefficient. Similarly, if $X_i \in {{\mathscr {S}}}$ is a product of ${\textbf{R}}_j \in {{\mathscr {R}}}$, the $X_i$ is a product of $\tilde{\textbf{R}}_j \in \tilde{{\mathscr {R}}}$ with the same stoichiometric coefficient.

We next provide a formal definition of the elementary modifications in Table 1.

Table 1 A list of elementary graph modifications studied in this paper

Definition 3

Let ${\mathscr {N}}=({\mathscr {S}},{\mathscr {R}})$ be a given BIN. We say that $\tilde{{\mathscr {N}}}:=(\tilde{{\mathscr {S}}},\tilde{{\mathscr {R}}})$ is an elementary modification of ${\mathscr {N}}$ if it satisfies one of the following statements:

• 1. (Reversal of a reaction) $\tilde{{\mathscr {S}}}= {{\mathscr {S}}}$, and $\exists {\textbf{R}}_j \in {{\mathscr {R}}}$ such that $\tilde{{\mathscr {R}}}={{\mathscr {R}}}\cup \{{\textbf{R}}_{-j}\}$.

• 2. (Adding an intermediate) $\tilde{{\mathscr {S}}}={{\mathscr {S}}}\cup \{X^*\}$, and $\exists {\textbf{R}}_j \in {{\mathscr {R}}}$ (written as ${\textbf{R}}_j=\sum _{i}\alpha _{ij}X_i \rightarrow \sum _i \beta _{ij} X_i$) such that $\tilde{{\mathscr {R}}}=({{\mathscr {R}}}/\{{\textbf{R}}_j\}) \cup \{\tilde{\textbf{R}}_{j},\tilde{\textbf{R}}^*\}$ where $\tilde{\textbf{R}}_j:= (\sum _{i}\alpha _{ij}X_i \rightarrow X^*)$, and $(\tilde{\textbf{R}}^*:= X^* \rightarrow \sum _i \beta _{ij} X_i)$.

• 3. (External Regulation) $\tilde{{\mathscr {S}}}= {{\mathscr {S}}}$, and $\exists X_k \in {{\mathscr {S}}}$ such that $\tilde{{\mathscr {R}}}={{\mathscr {R}}}\cup \{ X_k \leftrightharpoons \emptyset \}$.

• 4. (Conserved Regulation) $\tilde{{\mathscr {S}}}= {{\mathscr {S}}}\cup \{X_{n+1}\}$, and $\exists X_k \in {{\mathscr {S}}}$ such that $\tilde{{\mathscr {R}}}={{\mathscr {R}}}\cup \{ X_k \leftrightharpoons X_{n+1} \}$.

• 5. (Adding a feedback species) $\tilde{{\mathscr {S}}}= {{\mathscr {S}}}\cup \{X^*\}$, and $\exists {\textbf{R}}_j,{\textbf{R}}_k \in {{\mathscr {R}}}$ such that $\tilde{{\mathscr {R}}}=({{\mathscr {R}}}/\{{\textbf{R}}_j,{\textbf{R}}_k\}) \cup \{\tilde{\textbf{R}}_{j},\tilde{\textbf{R}}_k\}$ where $\tilde{\textbf{R}}_j$ is an extension of ${\textbf{R}}_j$ with $X^*$ as an extra product, and $\tilde{\textbf{R}}_k$ is an extension of ${\textbf{R}}_k$ with $X^*$ as an extra reactant.

• 6. (Adding a catalyst) $\exists X_i \in {{\mathscr {S}}}$ such that $\tilde{{\mathscr {S}}}= {{\mathscr {S}}}\cup \{X_i^-\}$, $\vert {\mathscr {R}}\vert =\vert \tilde{{\mathscr {R}}}\vert$, and every reaction $\tilde{\textbf{R}}_j\in \tilde{{\mathscr {R}}}$ is an extension of a corresponding reaction ${\textbf{R}}_j \in {{\mathscr {R}}}$. Furthermore, $X_i^-$ is a product of a reaction ${\textbf{R}}_j$ iff $X_i$ is a reactant of ${\textbf{R}}_j$ with the same stoichiometry coefficient, and $X_i^-$ is a reactant of a reaction ${\textbf{R}}_j$ iff $X_i$ is a product of ${\textbf{R}}_j$ with the same stoichiometry coefficient.

• 7. (Adding a dimer) $\exists X_i \in {{\mathscr {S}}}$ such that $\tilde{{\mathscr {S}}}= {{\mathscr {S}}}\cup \{X_i^+\}$, $\vert {\mathscr {R}}\vert =\vert \tilde{{\mathscr {R}}}\vert$, and every reaction $\tilde{\textbf{R}}_j\in \tilde{{\mathscr {R}}}$ is an extension of a corresponding reaction ${\textbf{R}}_j \in {{\mathscr {R}}}$. Furthermore, $\exists X_i \in {{\mathscr {S}}}$ such that $X_i^+$ is a reactant of a reaction ${\textbf{R}}_j$ iff $X_i^+$ is a reactant of ${\textbf{R}}_j$ with the same stoichiometry coefficient, and $X_i^+$ is a product of a reaction ${\textbf{R}}_j$ iff $X_i$ is a product of ${\textbf{R}}_j$ with the same stoichiometry coefficient.

Finally, a network $\tilde{{\mathscr {N}}}$ is a modification of ${{\mathscr {N}}}$ if it is a result of several elementary modifications. More formally:

Definition 4

A network $\tilde{{\mathscr {N}}}$ is a modification of ${{\mathscr {N}}}$ if there exists a finite sequence of networks ${{\mathscr {N}}}_0,{{\mathscr {N}}}_1,..,{{\mathscr {N}}}_q$, with ${{\mathscr {N}}}_0:={{\mathscr {N}}}, {{\mathscr {N}}}_{q}:=\tilde{{\mathscr {N}}}$, and for each $i\in {1,..,q}$, ${{\mathscr {N}}}_{i}$ is an elementary modification of ${{\mathscr {N}}}_{i-1}$.

In the subsequent sections, we provide results on modifications that preserve the stability of a given BIN.

Remark 2

The standard enzymatic catalysis reaction is a combination of three elementary modifications which are adding an intermediate, reversal, and then adding a catalyst. In other words, the reaction $S \rightarrow P$ is modified into $S \rightarrow C \rightarrow P$, then to $S \rightleftharpoons C \rightarrow P$, then to $S+E \leftrightharpoons C \rightarrow P+E$.

4.2 Linear networks with a sum-of-currents RLF

It is easy to see that the first few modifications in Table 1 are stability preserving when applied to a linear BIN. This is stated below.

Theorem 3

Let ${{\mathscr {N}}}$ be a given linear BIN, and let $\tilde{{\mathscr {N}}}$ be its modification generated by a finite sequence of elementary modifications that are limited to reversal of a reaction, adding an intermediate, external regulation of a species, and conserved regulation of a species. Then, V (5) is an RLF for $\tilde{{\mathscr {N}}}$.

Proof

The resulting network $\tilde{{\mathscr {N}}}$ after the application of the elementary modifications mentioned in the statement of the theorem is linear. Hence, the statement follows by Theorem 1. $\square$

The last two modifications in Table 1 are more interesting since they can modify a linear network into a nonlinear one. Nevertheless, we show that the resulting modified BIN continues to have an SoC RLF. The proof is provided in the appendix.

Theorem 4

Let ${{\mathscr {N}}}=({{\mathscr {S}}},{{\mathscr {R}}})$ be a given linear BIN, and let $\tilde{{\mathscr {N}}}=(\tilde{{\mathscr {S}}},\tilde{{\mathscr {R}}})$ be its modification generated by a finite sequence of elementary modifications that are limited to adding a catalyst and adding a dimer. Then, $V=\sum _{i=1}^{\vert {\mathscr {S}}\vert }\vert \dot{x}_i\vert$ is an RLF for $\tilde{{\mathscr {N}}}$.

Several modifications can be combined to yield enzymatic catalysis reactions (see Remark 2). Therefore, we can state the following corollary:

Corollary 5

Let ${{\mathscr {N}}}$ be a given linear BIN, and let $\tilde{{\mathscr {N}}}$ be its modification generated replacing linear reactions of the form $X_i \rightarrow X_j$, by nonlinear reactions of the form $X_i+E_{ij} \rightleftharpoons C_{ij} \rightarrow X_j + E_{ij}$. Let ${{\mathscr {S}}}_2$ be set of all the extra intermediates written as $C_{ij}$. Then, the function $V=\sum _{i=1}^{\vert {{\mathscr {S}}}\vert } \vert \dot{x}_i\vert + \sum _{C_{ij}\in {{\mathscr {S}}}_2} \vert \dot{c}_{ij}\vert$ is an RLF for $\tilde{{\mathscr {N}}}$.

Proof

The proof follows by using Theorem 3 for adding an intermediate and then reversal, i.e., modifying $X_i\rightarrow X_j$ to $X_i \leftrightharpoons C \rightarrow X_j$. Then, Theorem 4 to get the reaction $X_i +E_{ij}\leftrightharpoons C_{ij} \rightarrow X_j+E_{ij}$. $\square$

4.3 Networks with a max–min RLF

Networks that have a Max–Min RLF admit a different set of stability-preserving modifications as we show next. Note that the original BIN does not need to be linear as is stated in the following result.

Theorem 6

Let ${{\mathscr {N}}}$ be a BIN that admits a Max–Min RLF, and let $\tilde{{\mathscr {N}}}$ be its modification generated by a finite sequence of elementary modifications that are limited to adding an intermediate, adding a feedback species, adding a dimer, and adding a catalyst. Then, (6) is an RLF for $\tilde{{\mathscr {N}}}$.

Proof

Using the characterization in Theorem 2, any combination of the graph modifications mentioned in the statement of theorem do not create new independent vectors in the kernel of the stoichiometry matrix (i.e., it does not create new fluxes), and they do not make a single species a reactant in multiple reactions. Therefore, Theorem 2 applies to $\tilde{{\mathscr {N}}}$. $\square$

We study next the case of reversal. Since our formalism treats a reversible reaction as two reactions ${\textbf{R}}_j,{\textbf{R}}_{-j}$, then reversal of a reaction increases the number of fluxes, and hence violates the conditions required by Theorem 2. Nevertheless, as shown in [11], the result can be extended. We state the result here in the language of graph modifications:

Theorem 7

([11]) Let ${{\mathscr {N}}}=({{\mathscr {S}}},{{\mathscr {R}}})$ be a given network that satisfies the conditions of Theorem 2. Let ${{\mathscr {R}}}_r \subset {{\mathscr {R}}}$ be defined as follows: ${\textbf{R}}^* \in {{\mathscr {R}}}_r$ iff for each $X_i\in {{\mathscr {S}}}$ that is a product of ${\textbf{R}}^*$, $X_i$ is not a product of another reaction. Then, let $\tilde{{\mathscr {N}}}$ be a modification of ${{\mathscr {N}}}$ generated by the reversal of the reactions in ${{\mathscr {R}}}_r$. Then, (6) is an RLF for $\tilde{{\mathscr {N}}}=(\tilde{{\mathscr {S}}},\tilde{{\mathscr {R}}})$ where ${\mathcal {R}}(x)=\{R_{j}(x)-R_{-j}(x)\vert j=1,..,\vert {{\mathscr {R}}}\vert \}$, where $R_{-j}:\equiv 0$ if ${\textbf{R}}_{-j}\not \in \tilde{{\mathscr {R}}}$.

In addition, we can strengthen Corollary 5 to include modifications by processive enzymatic cycles [41]:

Corollary 8

Let ${{\mathscr {N}}}$ be a given BIN satisfying the conditions of Theorem 2, and let $\tilde{{\mathscr {N}}}$ be its modification generated by replacing reactions of the form $\sum _{i} \alpha _{ij} X_i \rightarrow \sum _{i} \beta _{ij} X_i$, by reactions of the form $\sum _{i} \alpha _{ij} X_i+E^* \rightleftharpoons C_{0}^* \rightleftharpoons C_{1}^* \rightleftharpoons .... \rightleftharpoons C_{m}^* \rightarrow \sum _{i} \beta _{ij} X_i + E^*$ for some positive integer m. Then, $\tilde{{\mathscr {N}}}$ admits a Max–Min RLF.

Proof

The statement can be proven by applying enzymatic catalysis (as in Remark 2) to get $\sum _{i} \alpha _{ij} X_i+E^* \rightleftharpoons C_{0}^* \rightarrow \beta _{ij} X_i + E^*$, then by the addition of intermediates $C_1^*,.. C_m^*$ and then reversals to get the required reaction. $\square$

5 Global stability and robust non-degeneracy

5.1 Global stability

Since our RLFs are non-strict, we need to verify global stability. A popular way is via LaSalle’s invariance principle. In our setting, a network ${{\mathscr {N}}}$ that admits an RLF V is said to satisfy the LaSalle’s principle if the following implication always holds: If a bounded solution $\tilde{x}(t)$ of (2) satisfies $\tilde{x}(t) \in \ker \dot{V}$ for all $t\ge 0$, then $\tilde{x}(t) \in \ker V$ for all $t \ge 0$, i.e., $V(\tilde{x}(t)))=0$.

Theorem 9

([11, 12]) Let a network ${{\mathscr {N}}}$ be given. Assume that ${{\mathscr {N}}}$ admits an RLF and it satisfies the LaSalle’s principle. Then,

• 1. Each bounded trajectory converges to the set of steady states,

• 2. if all the trajectories are bounded and there exists an isolated steady state relative to its stoichiometric class, then it is globally asymptotically stable.

5.1.1 Networks that admit an SoC RLF

In [11], an iterative algorithm has been proposed to check LaSalle’s invariance principle. However, in the next result, we show that it always holds for networks that satisfy Theorems 3 or 4.

Theorem 10

Let ${\mathscr {N}}$ be a linear network or a modification of a linear network that satisfies the conditions of Theorems 3 or 4. Then, it satisfies the LaSalle’s principle.

For linear networks, the statement has been shown in [39]. It remains to prove that generalization of the result to any nonlinear network that is a modification of a linear network. The proof is included in the Appendix.

5.1.2 Networks that admit a max–min RLF

Verification of LaSalle’s invariance principle for networks that admit Max–Min RLFs has been provided in [11] via a simple graphical condition. In order to introduce the next result, we need a definition. Consider a network ${{\mathscr {N}}}=({{\mathscr {S}}},{{\mathscr {R}}})$, then a reaction ${\textbf{R}}_k \in {{\mathscr {R}}}$ is said to be an ancestor of ${\textbf{R}}_j$ if there is a direct path from ${\textbf{R}}_k$ to ${\textbf{R}}_j$ on the Petri-net corresponding to ${{\mathscr {N}}}$. The result is stated in the following theorem:

Theorem 11

([11]) Let ${\mathscr {N}}$ be a network that satisfies the conditions of Theorems 2, 6, or 7. Then, ${{\mathscr {N}}}$ satisfies the LaSalle’s principle: if ${{\mathscr {N}}}$ is conservative, or if every pair of reactions share an ancestor.

5.2 Robust non-degeneracy

5.2.1 Definitions and review

In the previous subsection, we have shown that the trajectories converge to the set of steady states. However, existence of a steady state in a stoichiometric class does not automatically imply that it is isolated. Therefore, we study here the robust non-degeneracy of the Jacobian of (2) which can be written as $\Lambda :=\Gamma \partial R/\partial x =\Gamma V$, where $V\in {\mathcal {K}}_{{\mathscr {N}}}$. However, as mentioned in Sect. 2, the presence of a conservation law means that the positive orthant is a foliage of forward invariant stoichiometric classes. Therefore, the relevant entity for analysis is the reduced Jacobian $\Lambda _r$ which can be defined as follows. For a given $\Gamma \in {\mathbb {R}}^{n\times \nu }, V \in {\mathbb {R}}^{\nu \times n}$, denote $r:=\text{ rank }(\Gamma )$. Let $\{d_1,..,d_{n-r}\}$ be linearly independent left null vectors of $\Gamma$. In order to get a basis of ${\mathbb {R}}^n$, we add vectors to get the basis: $\{d_1,..d_{n-r},d_{n-r+1},..,d_{n}\}$, and get the transformation matrix:

$$\begin{aligned} T = [ d_1,...,d_n]^T. \end{aligned}$$

The Jacobian $\Lambda =\Gamma V$ in the new coordinates can be written as follows:

$$\begin{aligned} T \Gamma V T^{-1} = \begin{bmatrix} \Lambda _r &{} \Lambda _2 \\ 0 &{} 0 \end{bmatrix} \end{aligned}$$

(7)

The matrix $\Lambda _r \in {\mathbb {R}}^{r \times r}$ is the reduced Jacobian, and it is the Jacobian for the dynamics restricted to the stoichiometric class.

We are interested in its non-singularity for any admissible kinetics. Hence, we provide the following definition:

Definition 5

A network $({{\mathscr {S}}},{{\mathscr {R}}})$ is said to be robustly non-degenerate iff the reduced Jacobian $\Lambda _r$ defined in (7) is non-singular for all matrices $V \in {\mathcal {K}}_{{\mathscr {N}}}$.

Although (7) is written with a specific transformation matrix T, it is obvious to see that the non-singularity of the reduced Jacobian is independent of the specific choice of the matrix T.

In order to study the reduced Jacobian, we will use the concept of the essential determinant $\text{ det}_{{ess}}(\Lambda )$ which is defined as the sum of all $r\times r$ principal minors of $\Lambda$. The characterization can be stated as follows:

Lemma 12

([24]) Let $\Gamma \in {\mathbb {R}}^{n\times \nu }$ and $V \in {\mathbb {R}}^{\nu \times n}$ be given. Let $r:=\text{ rank }(\Gamma )$. The reduced Jacobian $\Lambda _r$ defined in (7) is non-singular iff $\text{ det}_{ess}(\Lambda )=\text{ det}_{ess}(\Gamma V)\ne 0$.

Hence, instead of explicitly computing the reduced Jacobian, our strategy will be to verify that the sum of $r\times r$ principal minors of the full Jacobian is nonzero for any $V \in {\mathcal {K}}_{{\mathscr {N}}}$. Our task is eased by the special properties of networks admitting a PWL RLF. In [11], we have proved that every principal minor of the Jacobian is, in fact, nonnegative. To state it more formally, we have the following definition:

Definition 6

A network ${{\mathscr {N}}}$ is said to be robustly $P_0$ if the Jacobian $-\Gamma V$ is $P_0$ for all $V\in {\mathcal {K}}_{{\mathscr {N}}}$, i.e., all its principal minors are nonnegative.

Hence, the result can be stated as follows.

Lemma 13

([11]) Let ${{\mathscr {N}}}$ be a given network. If it admits a PWL RLF, then it is robustly $P_0$.

Therefore, using the last two lemmas, we immediately get the following corollary:

Corollary 14

Let ${{\mathscr {N}}}$ be a given network that is robustly $P_0$. Then, ${{\mathscr {N}}}$ is robustly non-degenerate iff for every $V\in {\mathcal {K}}_N$, there exists a positive $r\times r$ principal minor, where $r={\text{ rank }}(\Gamma )$.

5.2.2 Computational testing of robust non-degeneracy

It is possible to computationally check robust non-degeneracy by testing the Jacobian at a finite number of points [12, 14, 42]. In fact, we have shown that one point is sufficient:

Theorem 15

([12, 42]) Let ${{\mathscr {N}}}$ be a network that admits a PWL RLF, and let $\Gamma \in {\mathbb {R}}^{n\times \nu }$ be the stoichiometry matrix with rank r. If $\exists V^* \in {\mathcal {K}}_{{{\mathscr {N}}}}$ such that $-\Gamma V^*$ has a positive essential determinant, then $-\Gamma V^*$ has a positive essential determinant for all $V\in {\mathcal {K}}_{{\mathscr {N}}}$, i.e., ${{\mathscr {N}}}$ is robustly non-degenerate.

In the next subsection, we provide our main result in this subsection, which is a graphical method to verify robust non-degeneracy.

Remark 3

Although Theorem 15 is stated in [12, 42] for networks admitting a PWL RLF, the proof holds for any robustly $P_0$ network.

5.2.3 Main result

Instead of directly verifying the non-degeneracy of a large network, we study it graphically. In other words, we consider the network as a modification of a simpler network. Therefore, we state our result which is proved in the appendix.

Theorem 16

Let ${{\mathscr {N}}}$ be a given BIN which is robustly $P_0$. Assume that ${{\mathscr {N}}}$ is robustly non-degenerate, and let $\tilde{{\mathscr {N}}}$ be its modification generated by a finite sequence of elementary modifications that are limited to reversal of a reaction, adding an intermediate, external regulation of a species, conserved regulation of a species, adding a catalyst, and adding a dimer. Then, if $\tilde{{\mathscr {N}}}$ is robustly $P_0$, it follows that $\tilde{{\mathscr {N}}}$ is robustly non-degenerate.

Remark 4

Any network that admits a PWL RLF is robustly $P_0$ using Lemma 13. Hence, Theorem 16 can be coupled with Theorems 3, 4, 6, 7 to show robust non-degeneracy of the modified networks.

5.3 Review of the consequences of robust non-degeneracy

Robust non-degeneracy of the Jacobian gives us a quick way to verify several key properties of BINs. For completeness, we review them here.

Uniqueness of steady states Lemma 13 implies that any network that admits a PWL RLF has a $P_0$ Jacobian, which excludes multiple non-degenerate steady states in the same stoichiometric class [24, 43]. Hence, we get the following:

Theorem 17

([12, 42]) Consider a network ${{\mathscr {N}}}$ that admits a PWL RLF and is robustly non-degenerate. Then, every positive steady state is unique relative to its stoichiometric class.

Exponential stability The following result follows from the properties of PWL Lyapunov functions:

Theorem 18

([12, 14, 42]) Let ${{\mathscr {N}}}$ be a network that admits a PWL RLF and robustly non-degenerate, then every positive steady state is exponentially asymptotically stable.

Global stability Using previous results, it can be readily seen that robust non-degeneracy coupled with the LaSalle’s principle implies global stability. However, it has been shown [15] that this can be strengthened to the following:

Theorem 19

([15]) Suppose that a network ${{\mathscr {N}}}$ admits a PWL RLF and is robustly non-degenerate, then every positive steady state is globally asymptotically stable relative to its stoichiometric class.

Remark 5

The statements in this subsection assume the existence of a positive steady state. One way to exclude the existence of steady states on the boundary is via verifying persistence. In other words, we need to guarantee that all the trajectories that start from the positive orthant do not asymptotically approach its boundary. Graphical conditions for persistence have already been developed in [26] and they are easily applicable as we will see in the next section.

6 Persistence

6.1 Definitions and review of previous results

For systems that evolve on the positive orthant, persistence simply means non-extinction [44, 45]. In other words, if a trajectory starts in the interior of the positive orthant, then it will not approach the boundary asymptotically. More precisely, a trajectory $\varphi (t;x_\circ )$ of (2) is said to be persistent if it satisfies $\liminf _{t\rightarrow \infty }\varphi (t;x_\circ )\gg 0$ whenever $x_\circ \gg 0$. A BIN network ${{\mathscr {N}}}$ is said to be robustly persistent if the previous statement holds for all bounded trajectories and for all admissible kinetics. A graphical notion of robust persistence for BINs has been introduced in [26, 46] using the concept of siphon which we define next.

Definition 7

Let ${{\mathscr {N}}}=({{\mathscr {S}}},{{\mathscr {R}}})$ be a given network. Then, a non-empty set $P \subset {{\mathscr {S}}}$ is called a siphon iff each input reaction associated to a species in P is also an output reaction associated to a (possibly-different) species in P. A siphon is said to be trivial if it contains the support of a conservation law, and it is said to be critical otherwise.

The main result is as follows:

Theorem 20

([26, 46, 47]) Let ${{\mathscr {N}}}$ be a network that lacks critical siphons. Then, ${{\mathscr {N}}}$ is robustly persistent.

This motivates the following definition:

Definition 8

A network ${{\mathscr {N}}}$ that lacks critical siphon is said to be graphically persistent.

6.2 Main result

We show here that graphical persistence is conserved under many types of modifications. We start with a general result which proved in the appendix:

Theorem 21

Let ${\mathscr {N}}$ be a given BIN which is graphically persistent, and let $\tilde{{\mathscr {N}}}$ be its modification generated by a finite sequence of elementary modifications that are limited to reversal of a reaction, external regulation of a species, conserved regulation of a species, adding an intermediate, and adding a dimer. Then, $\tilde{{\mathscr {N}}}$ is also graphically persistent.

We next show that the last theorem can be expanded for the classes of networks studied in this paper:

Theorem 22

Let ${\mathscr {N}}$ be a given BIN which is graphically persistent.

• 1. Let ${{\mathscr {N}}}$ be linear, and let $\tilde{{\mathscr {N}}}$ be its modification using any of the modifications listed in Theorems 3 or 4, then ${{\mathscr {N}}}$ is graphically persistent.

• 2. Assume that ${{\mathscr {N}}}$ is conservative and that it admits a Max–Min RLF, then it is graphically persistent.

Part 1 of Theorem 22 follows from Theorem 21 except for the case of adding a catalyst which is proved in the Appendix. Part 2 of Theorem 22 is proved in [11, Theorem 13].

7 Applications

7.1 Post-translational modification (PTM) cycles

The PTM cycle model is standard in systems biology [48]. The long-term dynamics of the PTM cycle have been a subject of extensive study using several methods. This includes monotonicity [9, 18], and RLFs [11, 12, 14]. In this paper, we show that the stability properties of the PTM cycle can be interpreted graphically in terms of the basic reversible reaction:

$$\begin{aligned} S \rightleftharpoons P, \end{aligned}$$

(8)

where S denotes the substrate, and P denotes the product. This simple motif admits both an SoC RLF and a Max–Min RLF. Furthermore, it is conservative, robustly non-degenerate, satisfies the LaSalle’s condition. In addition, it lacks critical siphons; hence, it is persistent [26]. Therefore, it satisfies the following statement which we call $(\star )$:

$(\star )$ “Each proper stoichiometric class contains a unique globally exponentially stable positive steady state.”

We show next that these properties are inherited by the modifications of the simple reversible reaction above.

7.1.1 The single PTM

We consider the single PTM cycle:

$$\begin{aligned} S + E \leftrightharpoons C \mathop {\longrightarrow }\limits P+E, P + F \leftrightharpoons D \mathop {\longrightarrow }\limits S+E. \end{aligned}$$

(9)

As noted in Remark 2, the reaction $S \rightarrow P$ can be modified into an enzymatic catalysis reaction. Using Corollary 5, we get that the PTM cycle above admits an SoC RLF. Furthermore, using Corollary 8 we get that it also admits a Max–Min RLF. Theorem 22 implies that it lacks critical siphons. Hence, using the results in §5 it satisfies the statement $(\star )$.

7.1.2 The PTM star

We can consider other modifications to (8). By adding a finite number of conserving regulations on S, we get the following network which we call the linear star (depicted in Fig. 2):

$$\begin{aligned} S \rightleftharpoons P_1, S \rightleftharpoons P_2,..., S \rightleftharpoons P_n. \end{aligned}$$

(10)

Then, using Corollary 5, we get that the PTM star (3)–(4) (depicted in Fig. 1) admits an SoC RLF since it is formed by enzymatic catalysis modifications. In addition, Theorem 22 implies that it lacks critical siphons. Hence, it satisfies that statement $(\star )$. Furthermore, to answer the question posed in the introduction. We can add the external regulation $\emptyset \leftrightharpoons S$ to (10), and then apply enzymatic catalysis to all other reactions to certify the existence of an SoC RLF. Since the network is no longer conservative, it satisfies the following statement: If a proper stoichiometric class contains a steady state, then it is a unique globally exponentially stable positive steady state.

Fig. 2

Stability of the linear star (left) implies stability of the PTM star (right). Using Corollary 5, existence of an RLF for the linear star implies the existence of an RLF for the PTM star depicted in Fig. 1. Details are provided in Sect. 7.1.2.

7.1.3 The processive multi-PTM cycle

Modifying (8) by adding intermediates gives the following network which we call the linear cycle (depicted in Fig. a):

$$\begin{aligned} S_0 \rightarrow S_1 \rightarrow ... \rightarrow S_n \rightarrow S_0, \end{aligned}$$

where $S_0:=S, S_n:=P$. Theorem 2 guarantees that the modified network has a Max–Min RLF. Corollary 8 implies that the following network admits a Max–Min RLF:

$$\begin{aligned} S_{i-1}+E_i&\rightleftharpoons C_{i1} \rightleftharpoons C_{i2} \rightleftharpoons .... \rightleftharpoons C_{im} \longrightarrow S_i + E_i, \\ S_{n}+E_n&\rightleftharpoons C_{n1} \rightleftharpoons C_{n2} \rightleftharpoons .... \rightleftharpoons C_{nm} \longrightarrow S_1 + E_n, \end{aligned}$$

$i=1,..,n-1$. The above network has been called the “all-encompassing” processive cycle, and its stability has been studied in [19] using monotone system techniques. Using our method, we show that the existence of an RLF follows by modifying the linear cycle (Fig. 3a) using processive enzymatic reactions to get the network depicted in Fig. 3b. In addition, Theorem 22 implies that it lacks critical siphons. Therefore, using the results in Sect. 5, it satisfies the statement $(\star )$.

7.1.4 The PTM chain

Consider now modifying (8) by a finite number of intermediates and reversals, we get the following network:

$$\begin{aligned} S_0 \rightleftharpoons S_1 \rightleftharpoons S_2... \rightleftharpoons S_n, \end{aligned}$$

(11)

where $S_0:=S, S_n:=P$. Corollary 5 implies that the following PTM chain admits an SoC RLF:

$$\begin{aligned} S_{i-1}+E_i&\rightleftharpoons C_{i} \longrightarrow S_i + E_i, \end{aligned}$$

(12)

$$\begin{aligned} S_{i-1}+F_i&\rightleftharpoons D_i \longrightarrow S_i + F_i, i=1,..,n. \end{aligned}$$

(13)

The existence of an SoC RLF of the PTM chain can be shown computationally for each given n by linear programming [12]. Nevertheless, Fig. 3c, d shows that the existence of an SoC RLF for each n follows from modifying a linear chain via enzymatic catalysis reactions. In addition, Theorem 22 implies that it lacks critical siphons. Therefore, using the results in §5, it satisfies the statement $(\star )$.

Fig. 3

Constructing an RLF for a nonlinear network from a linear one. a The linear cycle. b The processive multi-PTM cycle. The existence of an RLF follows from the existence of one for the linear cycle using Corollary 8. c The linear chain. d The PTM chain. The existence of an RLF follows from the existence of one for the linear chain using Corollary 5

Fig. 4

Additional examples for graphical RLF construction. a A linear BIN. b The McKeithan network. The existence of an RLF follows from the existence of one for the linear BIN in panel a using Corollary 5. c A one-directional linear chain. d The RFM. The existence of an RLF follows from the existence of one for the unidirectional linear chain using Corollary 5

7.2 T-cell kinetic proofreading

McKeithan [49] proposed a nonlinear BIN to explain T-cell’s ability to distinguish between different types of ligands. It is given as follows:

$$\begin{aligned} R+L&\rightleftharpoons C_0 \rightarrow C_1 \rightarrow ... \rightarrow C_n \nonumber \\ C_1&\rightarrow R+L, C_2 \rightarrow R+L,..., C_n \rightarrow R+L. \end{aligned}$$

(14)

Sontag [8] has studied the stability of the network using the theory of complex balance, while we have studied the network using computational RLF construction [12]. Here, we show that a stability certificate can be constructed by considering the network as a modification of a linear network. By noting that the species L is a dimer in the language of Table 1, we can see that (15) is a modification of the following network by the addition of a dimer:

$$\begin{aligned} RL&\rightleftharpoons C_0 \rightarrow C_1 \rightarrow ... \rightarrow C_n \nonumber \\ C_1&\rightarrow RL, C_2 \rightarrow RL,..., C_n \rightarrow RL. \end{aligned}$$

(15)

Hence, existence of an SoC RLF for (15) follows from Theorem 4. Figure a shows the linear network, while Fig. 4b shows the corresponding modified nonlinear network.

The set of steady states is globally stable by Theorem 10. We can also show robust non-degeneracy graphically as follows. We consider first a linear cycle $RL \rightarrow C_0 \rightarrow .. \rightarrow C_n \rightarrow RL$ which is robustly non-degenerate since it is a modification of $RL \rightleftharpoons C_n$. Then, adding reactions of the form $C_i \rightarrow RL$ would not increase the rank of the stoichiometry matrix; hence, the network in (15) is robustly non-degenerate using the same argument used in the proof of item 1 in Theorem 16. Finally, (14) is a modification of (15) by the addition of a dimer. Hence, robust non-degeneracy of (14) follows from Theorem 16. In addition, Theorem 22 implies that it lacks critical siphons. Therefore, using the results in §5, it satisfies the statement $(\star )$ for any N.

7.3 The ribosome flow model

The ribosome flow model (RFM) is a nonlinear system model of the process of translation initiation and elongation where it describes ribosome binding to codons on an mRNA that is being translated [50]. It has been shown [12] that the corresponding ODE can be written as a BIN with species $X_i, Y_i$ where $X_i$ is occupancy of the ith codon, while $Y_i$ is the vacancy of the ith codon. Hence, we get the following BIN (depicted in Fig. 4d):

$$\begin{aligned} Y_1&\rightarrow X_1, X_n \rightarrow Y_n, \\ X_i+Y_{i+1}&\rightarrow X_{i+1}+Y_{i}, \, i=1,..,n-1. \end{aligned}$$

The stability of the above network has been studied via monotonicity methods [20]. For a given n, the existence of an SoC RLF can be verified via linear programming [12]. Nevertheless, Fig. 4c, d shows that an SoC RLF can be constructed by merely noticing that the RFM is a modification generated by adding catalysts to the following unidirectional linear chain network: (depicted in Fig. 4c)

$$\begin{aligned} \emptyset \rightarrow X_1 \rightarrow X_2 \rightarrow ... \rightarrow X_n \rightarrow \emptyset . \end{aligned}$$

The same graphical technique can be applied to RFMs interconnected via a pool [21] (as Fig.   shows), or via multiple pools [22]. In addition, Theorem 22 implies that they lack critical siphons. Therefore, using the results in §5, and by noting that they lack critical siphons, all the aforementioned RFM variants satisfy the statement $(\star )$.

Fig. 5

Graphical construction of an RLF for the RFM with a pool [21]. a The linear network. b The corresponding modified nonlinear network. Stability follows via Corollary 5

8 Conclusion

In this work, we have proposed a graphical method to certify the existence of an RLF for a given network by reducing it via a certain set of admissible modifications to a network that is known to admit an RLF. Furthermore, our method can directly show that the stability of a given network is preserved under certain graph modifications. In addition, we have shown that properties of the original network such as global stability, robust non-degeneracy, and graphical persistence are invariant under such modifications. Using our methods, complex nonlinear networks of arbitrary size and arbitrary number of nonlinear reactions can be reduced into tractable networks.

Appendix: Proofs

1.1 Proof of Theorem 1

The function $V(x)={{\tilde{V}}}(R(x))$ is piecewise linear in terms of the rates; therefore, there exists a positive integer m such that the space ${\mathbb {R}}_{\ge 0}^\nu$ can be partitioned into non-empty-interior regions $\{{{\mathcal {W}}}_k\}_{k=1}^m \subset {\mathbb {R}}_{\ge 0}^\nu$ for which ${{\tilde{V}}}$ is linear on each of them and each region corresponds to a specific sign pattern for $\dot{x}$. The geometry of such partition is discussed more thoroughly in [11].

Fix k. There exists $c_{ij}^{(k)}, c_i^{(k)}, \delta ^{(k)}$ such that:

$$\begin{aligned} V(x)&=\sum _{\begin{array}{c} i,j\\ i\ne j \end{array}} c_{ij}^{(k)} R_{ij}(x_i) + \sum _{i} c_{i}^{(k)} R_{i}(x_i) + c_0^{(k)} \nonumber \\&=:c^{(k)^T}R(x), \ R(x) \in {{\mathcal {W}}}_k. \end{aligned}$$

(16)

Since V is defined as the $\ell _1$ norm of $\dot{x}$, then the sign of $\dot{x}$ is constant and nonzero on ${{\mathcal {W}}}_k^\circ$. Therefore, we denote $\sigma _i:={\text{ sgn }}(\dot{x}_i) \in \{\pm 1\}$ on ${{\mathcal {W}}}_k^\circ$, where the superscript “$\circ$” denotes the interior of a set.

We claim that each term in the expression (16) has a nonpositive Lie derivative on ${{\mathcal {W}}}_k^\circ$. In order to show that, we first examine terms of the form $c_{ij}^{(k)} R_{ij}(x_i)$ where $c_{ij}^{(k)} \ne 0$ for some i, j. We will show that $c_{ij}^{(k)} \dot{R}_{ij}(x_i) \le 0$ for $R(x) \in {{\mathcal {W}}}_k^\circ$. As evident by examining (5), the reaction rate $R_{ij}$ appears only in $\dot{x}_i$ with coefficient $-1$ and in $\dot{x}_j$ with coefficient $+1$. W.l.o.g, assume that $c_{ij}^{(k)}>0$. There are four possible combinations $\sigma _i^{(k)},\sigma _j^{(k)}>0$, $\sigma _i^{(k)},\sigma _j^{(k)}<0$, $\sigma _i^{(k)}>0,\sigma _j^{(k)}<0$, and $\sigma _i^{(k)}<0,\sigma _j^{(k)}>0$. The first two give $c_{ij}^{(k)}=0$ and the third gives $c_{ij}^{(k)}=-2<0$. Hence, we conclude that $\sigma _i^{(k)}<0,\sigma _j^{(k)}>0$.

Therefore, ${\text{ sgn }}(c_{ij}^{(k)} \dot{R}_{ij}(x_i))\! =\! {\text{ sgn }}(c_{ij}^{(k)} (\partial R_{ij}(x_i)/\partial x_i) \dot{x}_i) \mathop {=}\limits ^{ } {\text{ sgn }}(c_{ij}^{(k)} \sigma _i^{(k)})\le 0$ for $R(x) \in {{\mathcal {W}}}_k^\circ$, where the last equality follows by the monotonicity of $R_{ij}$.

Next, we examine $c_i^{(k)} R_i(x_i)$ for some i where $c_i^{(k)} \ne 0$. W.l.o.g, assume that $c_i^{(k)}>0$. Since $R_i$ appears only in $\dot{x}_i$ with coefficient $-1$, then $\sigma _i^{(k)}<0$. Therefore, ${\text{ sgn }}(c_{i} \dot{R}_{i}(x_i)) = {\text{ sgn }}(c_{i}^{(k)} (\partial R_{i}(x_i)/\partial x_i) \dot{x}_i) \mathop {=}\limits ^{ } {\text{ sgn }}(c_{i}^{(k)} \sigma _i)\le 0$ for $R(x) \in {{\mathcal {W}}}_k^\circ$.

Since k, i, j have been chosen arbitrarily, we conclude that $\dot{V}(x) \le 0$ whenever $R(x) \in {{\mathcal {W}}}_k^\circ$ for some k. It remains to show that $\dot{V}(x) \le 0$ when $R(x) \in \partial {{\mathcal {W}}}_k$ for some k where “$\partial$” denotes the boundary of a set. To that end, similar to [11][Proof of Theorem 2], the Dini’s derivative can be written as $\dot{V}(x)=\max _{k \in K_{x(t)}} {c^{(k)}}^T \dot{R}(x) \le 0$ where $K_{x(t)}=\{k \vert R(x) \in {{\mathcal {W}}}_k\}$. $\square$

1.2 Proof of Theorem 4

Let $\Gamma$ be the stoichiometry matrix for $({{\mathscr {S}}},{{\mathscr {R}}})$. Since the modifications are limited to adding a catalysis or adding a dimer, then every reaction in $\tilde{{\mathscr {R}}}$ is an extension of a corresponding reaction in ${{\mathscr {R}}}$. Hence, we can write ${{\tilde{\Gamma }}}=[\Gamma ^T,\Gamma _2^T]^T$ as the stoichiometry matrix for $({{\tilde{{{\mathscr {S}}}}}},{{\tilde{{{\mathscr {R}}}}}})$. Let $\dot{x}=\Gamma R(x), \dot{{{\tilde{x}}}}={{\tilde{\Gamma }}} {{\tilde{R}}}({{\tilde{x}}})$ be the corresponding ODEs. Hence, we can write ${{\tilde{x}}}=[x^T,x_2^T]^T$, where $x_2$ corresponds to the concentrations of the species in ${{\tilde{{{\mathscr {S}}}}}}/{{\mathscr {S}}}$.

Note that all the species in ${{\tilde{{{\mathscr {S}}}}}}/{{\mathscr {S}}}$ are either catalysts or dimers. We include an additional assumption to simplify the notation: For each species $X_i \in {{\mathscr {S}}}$, we assume that there exists at most one corresponding catalyst species in ${{\tilde{{{\mathscr {S}}}}}}/{{\mathscr {S}}}$, and it is denoted by $X_i^-$. Similarly, we assume that there exists at most one corresponding dimer species, and the corresponding species is denoted as $X_i^+$. The corresponding concentrations are $x_i, x_i^-, x_i^+$. The proof can be generalized easily without the last assumption.

Our construction implies that $\dot{x}_i=\dot{x}_i^+ = -\dot{x}_i^-$. Hence, $V(x)=0$ iff $\dot{x}=0$. Therefore, V is positive-definite. We next show that it is non-increasing.

Similar to the proof of Theorem 1, we consider a region ${{\mathcal {W}}}_k$ for which V is linear and has a fixed sign pattern for $\dot{x}$. Fix k, There exists $c_{ij}^{(k)}, c_i^{(k)}, \delta ^{(k)}$ such that:

$$\begin{aligned} V(x)&=\sum _{\begin{array}{c} i,j\\ i\ne j \end{array}} c_{ij}^{(k)} R_{ij}(x_i,x_i^+,x_j^-) + \sum _{i} c_{i}^{(k)} R_{i}(x_i,x_i^+) + c_0^{(k)} \nonumber \\ {}&=:c^{(k)^T}R(x), \ R(x) \in {{\mathcal {W}}}_k. \end{aligned}$$

(17)

We claim that each term in the expression (16) has a nonpositive Lie derivative on ${{\mathcal {W}}}_k^\circ$. In order to show that, we first examine $c_{ij}^{(k)} R_{ij}(x_i)$ where $c_{ij}^{(k)} \ne 0$ for some i, j. We will show that $c_{ij}^{(k)} \dot{R}_{ij}(x_i) \le 0$ for $R(x) \in {{\mathcal {W}}}_k^\circ$. Since the candidate RLF sums only the species in ${{\mathscr {S}}}$, the reaction rate $R_{ij}$ appears only in $\dot{x}_i$ with coefficient $-1$ and in $\dot{x}_j$ with coefficient $+1$. W.l.o.g, assume $c_{ij}^{(k)}>0$. Similar to the proof of Theorem 1, we get that $\sigma _i^{(k)}<0,\sigma _j^{(k)}>0$. Since $\sigma _i^{+(k)} = \sigma _i^{(k)}$, and $\sigma _j^{-(k)} = -\sigma _j^{(k)}$ we get

$$\begin{aligned} {\text{ sgn }}(c_{ij}^{(k)} \dot{R}_{ij}(x_i,x_i^+,x_j^-))&= {\text{ sgn }}\left( c_{ij}^{(k)} \left( \frac{\partial R_{ij} }{\partial x_i} \dot{x}_i + \frac{\partial R_{ij} }{\partial x_i^+} \dot{x}_i^+ + \frac{\partial R_{ij} }{\partial x_j^-} \dot{x}_j^- \right) \right) \\ {}&\mathop {=}\limits ^{ } {\text{ sgn }}(c_{ij}^{(k)} (\sigma _i^{(k)}+\sigma _i^{+(k)}+\sigma _j^{-(k)} ))\le 0 \end{aligned}$$

for $R(x) \in {{\mathcal {W}}}_k^\circ$, where the last equality follows by the monotonicity of $R_{ij}$.

Since k, i, j have been chosen arbitrarily, we can use the same arguments used in the proof of Theorem 1 to conclude that $\dot{V}(x)\le 0$ for all x. $\square$

1.2.1 Proof of Theorem 10

Let $\tilde{{\mathscr {N}}}=(\tilde{{\mathscr {S}}},\tilde{{\mathscr {R}}})$ be a modification of a linear network ${{\mathscr {N}}}=({{\mathscr {S}}},{{\mathscr {R}}})$ by adding a catalysts or dimers. We use the standard LaSalle’s principle [35]. Let x(t) be a trajectory of (2) that is contained in $\ker \dot{V}$, i.e., $\dot{V}(x(t))\equiv 0$, and $V(x(t))\equiv {{\bar{V}}}_1 \ge 0$. In order to prove global stability, we need to show that ${{\bar{V}}}_1=0$.

Recall that $V(x)=\sum _{i=1}^{\vert {{{\mathscr {S}}}}\vert } \vert {\dot{x}_i}\vert$, and $\sigma _i(t)={\text{ sgn }}(x_i(t))$. Then, let the time-dependent sets $\Sigma _+,\Sigma _-,{{\bar{\Sigma }}}_+,{{\bar{\Sigma }}}_- \subset \{1,..,\vert {{{\mathscr {S}}}}\vert \}$ be defined as: $\Sigma _+(t)=\{i\vert \sigma _i(t)> 0\}$, $\Sigma _-(t)=\{i\vert \sigma _i(t) < 0\}$, ${{\bar{\Sigma }}}_+(t)=\{i\vert \sigma _i(t)\ge 0\}$, and ${{\bar{\Sigma }}}_-(t)=\{i\vert \sigma _i(t)\le 0\}$.

Since $V(x)=\sum _{i=1}^{\vert {{{\mathscr {S}}}}\vert }\vert {\dot{x}_i}\vert$, we can write:

$$\begin{aligned} V(x(t)) = \sum _{i\in \Sigma _+(t)} \dot{x}_i(t) - \sum _{i\in \Sigma _-(t)} \dot{x}_i(t) \equiv {{\bar{V}}}_1. \end{aligned}$$

(18)

If either one of the sets $\Sigma _+(t), \Sigma _-(t)$ is empty for some t, then this implies that $\dot{x}(t)\equiv 0$, and hence ${{\bar{V}}}_1=0$ which proves the statement. Therefore, we assume that both $\Sigma _+(t), \Sigma _-(t)$ are non-empty for all t.

Similar to the proof of Theorem 4, let $x_i^+$ denote the concentration of the dimer that corresponds to $X_i$, and let $x_j^-$ be the concentration of the dimer that corresponds to $X_j$. Denote $\dot{R}_{ij}(x_i,x_i^+,x_j^-):= \frac{\partial R_{ij}}{\partial x_i}\dot{x}_i + \frac{\partial R_{ij}}{\partial x_i^+}\dot{x}_i^+ + \frac{\partial R_{ij}}{\partial x_j^-}\dot{x}_j^-, \dot{R}_{i}(x_i,x_i^+):= \frac{\partial R_{i}}{\partial x_i}\dot{x}_i+\frac{\partial R_{i}}{\partial x_i^+}\dot{x}_i^+.$ Using the argument in the Theorem 4, a term of the form $R_{ij}(x)$ appears in V with a positive coefficient only if $\sigma _i\ge 0$, $\sigma _j\le 0$, and $\sigma _i \sigma _j\ne 0$. Similarly, $R_{ij}(x)$ appears in V with a negative coefficient only if $\sigma _i\le 0$, $\sigma _j\ge 0$, and $\sigma _i \sigma _j\ne 0$. Hence, we can write:

$$\begin{aligned} \dot{V}= \sum _{(i,j)\in {{\bar{\Sigma }}}_- \times {{\bar{\Sigma }}}_+} \left( \rho _{ij} \dot{R}_{ij}(x) + \dot{R}_i(x) \right) - \sum _{(i,j)\in {{\bar{\Sigma }}}_+ \times {{\bar{\Sigma }}}_-}\left( \rho _{ij} \dot{R}_{ij}(x) - \dot{R}_i(x) \right) \equiv 0, \end{aligned}$$

for some $\rho _{ij}>0$. Note that the dependence on t in the equation above has been dropped for notational brevity.

As in the proof of Theorem 4, each term is nonpositive. Hence, $\dot{V}\equiv 0$ implies that each term is identical to zero. We make several conclusions from the last statement:

First, $\forall i \in \Sigma _+ \cup \Sigma _-, \dot{R}_i(x(t))\equiv 0$. Furthermore, by definition, $\dot{R}_i(x(t))=0$ for $i \not \in \Sigma _+ \cup \Sigma _i$. Therefore, we get that $\forall i, \dot{R}_i(x(t))\equiv 0$. Hence, $\sum _{i=1}^{\vert {{{\mathscr {S}}}}\vert } \dot{x}_i(t) = \sum _{i=1}^{\vert {{{\mathscr {S}}}}\vert } \sum _{j\ne i} (\dot{R} _{ji}(x(t))-\dot{R}_{ij}(x(t))) \equiv 0.$ Therefore, we get $\sum _{i=1}^{\vert {{{\mathscr {S}}}}\vert } \dot{x}_i(t) \equiv {{\bar{V}}}_2$ for some constant ${{\bar{V}}}_2$. Hence,

$$\begin{aligned} V(x(t)) = \sum _{i\in \Sigma _+(t)} \dot{x}_i(t)+ \sum _{i\in \Sigma _-(t)} \dot{x}_i(t) \equiv {{\bar{V}}}_2. \end{aligned}$$

(19)

Second, fix $i \in \Sigma _+(t)$. Then, for all $j\in \Sigma _-(t)$, we have $\dot{R}_{ij}(x(t)) \equiv 0$. Hence, we claim the following: If $i \in \Sigma _i^+(t')$ for some $t'>0$, then $i \in \Sigma _i^+(t')$ for all $t>t'$. To show this, we can write

$$\begin{aligned} \ddot{x}_i(t) = -F_i(\dot{x}_i) + \overbrace{ \sum _{\begin{array}{c} i \in \Sigma _+ \\ j\ne i \end{array}} \dot{R}_{ji}(x)}^{>0} + \overbrace{\sum _{j\in \Sigma _-} {\dot{R}}_{ji}(x)}^{=0}, \end{aligned}$$

where $F_i$ lumps all the terms that corresponds to reactions for which $X_i$ is a reactant. It can be noted immediately that $F_i(\dot{x}_i)>0$ is positive when $\dot{x}_i>0$, and $F_i(0)=0$. The second term is positive since $R_{ji}$ is monotone and $\dot{x}_j>0$. The third term is zero as we have shown earlier. Therefore, $\left. \ddot{x}_i\right| _{\dot{x}_i=0} >0$. Hence, $\dot{x}_i(t)>0$ for all $t>t'$, i.e., $i \in \Sigma _+(t)$ for all $t\ge t'$ as claimed.

Similarly, if $i \in \Sigma _i^-(t')$ for some $t'>0$, then $i \in \Sigma _i^-(t')$ for all $t>t'$. Since $\Sigma _+,\Sigma _- \subset \{1,..,\vert {{{\mathscr {S}}}}\vert \}$, there exists $T>0$ and constant sets $\Sigma _+^*, \Sigma _-^*$ with $\Sigma _+(t)=\Sigma _+^*,\Sigma _-(t)=\Sigma _-^*$ for all $t\ge T$. Therefore, we can write (18), (19) as:

$$\begin{aligned} \sum _{i\in \Sigma _+^*} \dot{x}_i(t)- \sum _{i\in \Sigma _-^*} \dot{x}_i(t) \equiv {{\bar{V}}}_1, ~\sum _{i\in \Sigma _+^*} \dot{x}_i(t)+ \sum _{i\in \Sigma _-^*} \dot{x}_i(t) \equiv {{\bar{V}}}_2. \end{aligned}$$

Hence, $\sum _{i\in \Sigma _+^*} \dot{x}_i(t) \equiv \tfrac{1}{2} ({{\bar{V}}}_1+{{\bar{V}}}_2)$, and $\sum _{i\in \Sigma _-^*} \dot{x}_i(t) \equiv \tfrac{1}{2} ({{\bar{V}}}_2-{{\bar{V}}}_1)$. By integrating the last two equations, we note that we can only have ${{\bar{V}}}_1={{\bar{V}}}_2=0$ since x(t) has been assumed to be bounded. $\square$

1.3 Proof of Theorem 16

Before proving the result, we need few preliminary lemmas.

Lemma 23

Let $\Lambda =-\Gamma V$, where $\Gamma \in {\mathbb {R}}^{n \times \nu }$, $V \in {\mathbb {R}}^{\nu \times R}$. Let $r={\text{ rank }}(\Lambda )$. Then,

• 1. Let I be an arbitrary subset of $\{1,..,n\}$ with $\vert I\vert =r$. The corresponding principal minor can be written as:

  $$\begin{aligned} \text{ det}_I( \Lambda )=\text{ det}_I (-\Gamma V) = \sum _{\begin{array}{c} J \subset \{1,..,\nu \}\\ \vert J\vert =r \end{array}} \text{ det }(-\Gamma _{IJ})\text{ det }(V_{JI}), \end{aligned}$$

  (20)

  where $\Gamma _{IJ},V_{JI}$ denote the submatrices of $\Gamma ,V$ with the row and column indices specified in I and J, respectively.

• 2. The essential determinant of $\Lambda$ can be written as:

  $$\begin{aligned} \text{ det}_{ess}(\Lambda ) = \sum _{\begin{array}{c} I \subset \{1,..,n\}\\ \vert I\vert =r \end{array}}\sum _{\begin{array}{c} J \subset \{1,..,\nu \}\\ \vert J\vert =r \end{array}} \text{ det }(-\Gamma _{IJ})\text{ det }(V_{JI}). \end{aligned}$$

  (21)

Proof

It follows immediately from the Cauchy–Binet’s formula [51] and the definition of the essential determinant. $\square$

Before stating the next lemma, we need some notation. Recall that ${\mathcal {K}}_{{\mathscr {N}}}$ is the set of all possible Jacobian matrices $\partial R/\partial x$ defined on ${\mathbb {R}}_+^n$. Hence, any $V \in {\mathcal {K}}_{{\mathscr {N}}}$ can have s nonzero entries which is equal to the number of reactant-reaction pairs. Let s be the number of nonzero entries of V, we list them as $v_1,..,v_s > 0$. Next, we show that each term in the expansion (20) is nonnegative. In other words,

Lemma 24

Let ${{\mathscr {N}}}$ be a network that is robustly $P_0$ and non-degenerate. Let the stoichiometry matrix be $\Gamma \in {\mathbb {R}}^{n\times \nu }$, and let $r={\text{ rank }}(\Gamma )$. Then, $\forall I\subset \{1,..,n\}, \forall J\subset \{1,..,\nu \}$ with $\vert I\vert =\vert J\vert =r$ we have $\text{ det }(-\Gamma _{IJ})\text{ det }(V_{JI}) \ge 0$ for all $V \in {\mathcal {K}}_{{\mathscr {N}}}$.

Proof

As noted before, $v_1,..,v_s$ are the nonzero entries of $V \in {\mathcal {K}}_{{\mathscr {N}}}$. Hence, we can think of $\text{ det }(V_{JI})$ as a polynomial in $v_1,..,v_s$. For the sake of contradiction, assume that there exists $I^*,J^*$ with $\vert {I}\vert ^*=\vert {J}\vert ^*=r$ such that $\text{ det }(-\Gamma _{I^*J^*})\text{ det }(V_{J^*I^*})<0$ for some $V\in {\mathcal {K}}_{{\mathscr {N}}}$. Hence, $\text{ det }(V_{J^*I^*})$ must have a monomial term $m^*=\prod _\ell v_\ell$ such that $\text{ det }(-\Gamma _{I^*J^*}) m^* <0$.

The corresponding principal minor can be written using (20) as $\text{ det}_{I^*}(-\Gamma V)=\sum _{J, \vert {J}\vert =r} \text{ det }(-\Gamma _{I^*J})\text{ det }(V_{JI^*})$. All the entries in V that do not appear in $m^*$ can be set to be arbitrarily small. Since the determinant is homogeneous in the entries of the corresponding matrix, all the terms other than $\text{ det }(-\Gamma _{I^*J^*}) m^* <0$ will become arbitrarily small in the expansion (20) which makes the principal minor $\text{ det}_{I^*}(-\Gamma V)<0$ for some $V \in {\mathcal {K}}_{{\mathscr {N}}}$. However, ${{\mathscr {N}}}$ is robustly $P_0$ which means that all principal minors of $-\Gamma V$ are nonnegative for any $V \in {\mathcal {K}}_{{\mathscr {N}}}$ which is a contradiction. $\square$

We are ready now to prove the statement of the theorem. Let $\tilde{{\mathscr {N}}}=(\tilde{{\mathscr {S}}},\tilde{{\mathscr {R}}})$ be the modified network, and let $\tilde{\Gamma }$ be its stoichiometry matrix. Furthermore, since ${{\mathscr {N}}}$ is robustly non-degenerate, and using Corollary 14 and Lemma 24, $\exists I\subset \{1,..,n\}, \exists J\subset \{1,..,\nu \}$ with $\vert {I}\vert =r$,$\vert {J}\vert =r$ such that $\text{ det }(-\Gamma _{IJ} )\text{ det }(V_{JI})>0$.

The proof is divided by the graph modification under consideration.

• 1. Reversal: Recall $\tilde{{\mathscr {S}}}={{\mathscr {S}}}$ and $\tilde{{\mathscr {R}}}={{\mathscr {R}}}\cup \{{\textbf{R}}_{-j}\}$ for some j. Let $\Gamma _j$ denote the jth column of $\Gamma$. Then, $\tilde{\Gamma }=[\Gamma ,-\Gamma _j]$. Hence, ${\text{ rank }}(\tilde{\Gamma })={\text{ rank }}(\Gamma )=r$. Using Corollary 14, we need to show the existence of a positive $r\times r$ principal minor. We consider the minor corresponding to I for $\tilde{{\mathscr {N}}}$. It can be seen immediately that is also positive since addition of a reverse of a reaction can only add new nonnegative terms (by Lemma 24) to the expansion (20).

• 2. Adding an intermediate: Recall that the network ${{\mathscr {N}}}$ is modified by replacing a reaction of the form ${\textbf{R}}_j=\sum _{i}\alpha _{ij}X_i \rightarrow \sum _i \beta _{ij} X_i$ by two reactions $\tilde{\textbf{R}}_j:= (\sum _{i}\alpha _{ij}X_i \rightarrow X_{n+1})$, and $(\tilde{\textbf{R}}_{\nu +1}:= X_{n+1} \rightarrow \sum _i \beta _{ij} X_i)$. W.l.o.g, assume that ${\textbf{R}}_j$ is the last reaction, i.e., $j=\nu$. Let $\Gamma _\nu$ be the last column of $\Gamma$, hence we write $\Gamma =[{{\hat{\Gamma }}},\Gamma _\nu ]$. Define the two vectors $\Gamma _\nu ^+,\Gamma _\nu ^-$ entry-wise as $[\Gamma _\nu ^+]_i=\max \{[\Gamma _\nu ]_i,0\}, [\Gamma _\nu ^-]_i=\min \{[\Gamma _\nu ]_i,0\}$, respectively. In other words, $\Gamma _\nu ^+$ contains the positive entries of $\Gamma _\nu ^+$ while $\Gamma _\nu ^-$ contains the negative entries of $\Gamma _\nu$. Hence, by construction, $\Gamma _{\nu }=\Gamma _\nu ^++\Gamma _\nu ^-$ Then, it can be seen that:

  $$\begin{aligned} \tilde{\Gamma }=\begin{bmatrix}{{\hat{\Gamma }}} &{} \tilde{\Gamma }_\nu ^- &{} \Gamma _\nu ^+ \\ 0 &{} 1 &{} -1 \end{bmatrix}. \end{aligned}$$

  Let $\tilde{I}:=I \cup \{n+1\}, \tilde{J}=J\cup \{\nu +1\}$. We will be computing $\text{ det }(-{{\tilde{\Gamma }}}_{\tilde{I}\tilde{J}})$ and $\text{ det }(\tilde{V}_{\tilde{I}\tilde{J}})$. We start with the latter. The new reaction has only one reactant which is $X_{n+1}$. Hence, we can write:

  $$\begin{aligned} \text{ det }(\tilde{V}_{\tilde{J}\tilde{I}})=\text{ det }\left( \begin{bmatrix} V_{JI} &{} 0 \\ 0 &{} 1 \end{bmatrix} \right) = \text{ det }(V_{JI}). \end{aligned}$$

  (22)

  Next, we consider $\text{ det }(-{{\tilde{\Gamma }}}_{\tilde{I}\tilde{J}})$. We study two cases: $\nu \in J$ and $\nu \not \in J$. $\underline{\hbox {Case 1:} \nu \in \hbox {J:}}$ Let $J^*=J/\{\nu \}$. We get:

  $$\begin{aligned} \text{ det }(-\Gamma _{\tilde{I} \tilde{J}})&=\text{ det }\left( -\begin{bmatrix} {{\hat{\Gamma }}}_{IJ^*} &{} {{\tilde{\Gamma }}}_{\nu ,I}^- &{} {{\tilde{\Gamma }}}_{\nu ,I}^+ \\ 0 &{} 1 &{} -1 \end{bmatrix}\right) \mathop =^{(\heartsuit )} \text{ det }\left( -\begin{bmatrix} {{\hat{\Gamma }}}_{IJ^*} &{} {{\tilde{\Gamma }}}_{\nu ,I} &{} {{\tilde{\Gamma }}}_{\nu ,I}^+ \\ 0 &{} 0 &{} -1 \end{bmatrix}\right) \nonumber \\&= \text{ det }\left( -\begin{bmatrix} \Gamma _{IJ} &{} {{\tilde{\Gamma }}}_{\nu ,I}^+ \\ 0 &{} -1 \end{bmatrix}\right) = \text{ det }(-\Gamma _{IJ}), \end{aligned}$$

  (23)

  where $(\heartsuit )$ follows by the fact that the determinant is invariant under the addition of the last two columns. $\underline{\hbox {Case 2:} \nu \not \in \hbox {J:}}$ We get:

  $$\begin{aligned} \text{ det }(-\Gamma _{\tilde{I} \tilde{J}})&=\text{ det }\left( -\begin{bmatrix} \Gamma _{IJ} &{} {{\tilde{\Gamma }}}_{\nu ,I}^+ \\ 0 &{} -1 \end{bmatrix}\right) = \text{ det }(-\Gamma _{IJ}), \end{aligned}$$

  (24)

  Hence, using (22), (23), (24), we get that $\text{ det }(-{{\tilde{\Gamma }}}_{\tilde{I}\tilde{J}} {{\tilde{V}}}_{\tilde{J}\tilde{I}})=\text{ det }(- \Gamma _{ I J} V_{ J I}) > 0$. Finally, since $\tilde{{\mathscr {N}}}$ admits an RLF, then Lemma 24 and Corollary 14 imply that the $(r+1)\times (r+1)$ principal minor corresponding to $\tilde{I}$ for modified network $\tilde{{\mathscr {N}}}$ is positive, hence $\tilde{{\mathscr {N}}}$ is robustly non-degenerate.

• 3. External regulation of a species: $\tilde{{\mathscr {S}}}= {{\mathscr {S}}}$, and $\exists X_i \in {{\mathscr {S}}}$ such that $\tilde{{\mathscr {R}}}={{\mathscr {R}}}\cup \{ X_i \leftrightharpoons \emptyset \}$. W.l.o.g, assume that $i=n$. Then we study two cases: ${\text{ rank }}(\tilde{\Gamma })={\text{ rank }}\Gamma$ and ${\text{ rank }}(\tilde{\Gamma })={\text{ rank }}\Gamma +1$. In first case, we let $\tilde{I}=I, \tilde{J}=J$. Hence, $\text{ det }(-\tilde{\Gamma }_{\tilde{I} \tilde{J}})\text{ det }(\tilde{V}_{\tilde{J} \tilde{I}})=\text{ det }(- \Gamma _{I J})\text{ det }( V_{ JI})>0$. Therefore, using Lemma 24 and Corollary 14, $\tilde{{\mathscr {N}}}$ is robustly non-degenerate. We now study the case in which ${\text{ rank }}\tilde{\Gamma }=1+{\text{ rank }}\Gamma =1+r$. We will first claim that it must be possible to choose I such that $n\not \in I$. As a proof, consider the contrary. Then, this means that $\forall I\subset \{1,..,n\}$ that satisfies $\vert {I}\vert =r$ and $n\not \in I$, we have $\text{ det}_I(-\Gamma V)=0$. Using Lemma 12, this means that removing $\gamma _n^T$ (the nth row of $\Gamma$), i.e., removing $X_n$ from ${{\mathscr {N}}}$, will cause the rank of $\Gamma$ to drop from r to $r-1$. Hence, this means that $\gamma _n^T$ is linearly independent from the other rows of $\Gamma$. However, adding the reactions $\{ X_i \leftrightharpoons \emptyset \}$ will only modify the nth row in $\Gamma$. Since $\gamma _n$ is already independent of the other rows of $\Gamma$, the rank cannot increase, which is a contradiction. Therefore, when ${\text{ rank }}\tilde{\Gamma }=r+1$, we let I be chosen such that $n\not \in I$. Hence, let $\tilde{I}=I \cup \{n\}, \tilde{J}=J\cup \{\nu +1\}$, where $\tilde{\textbf{R}}_{\nu +1}:= X_n \rightarrow \emptyset$. Therefore, we can write,

  $$\begin{aligned} \text{ det }(-\tilde{\Gamma }_{\tilde{I} \tilde{J}})&=\text{ det }\left( -\begin{bmatrix} \Gamma _{IJ} &{} 0 \\ 0 &{} -1 \end{bmatrix}\right) = \text{ det }(-\Gamma _{IJ}), \end{aligned}$$

  (25)

  Similarly, $\text{ det }(\tilde{V}_{\tilde{J} \tilde{I}})=\text{ det }( V_{ J I})$ using the same argument as in (22). Therefore, $\text{ det }(-\tilde{\Gamma }_{\tilde{I} \tilde{J}})\text{ det }(\tilde{V}_{\tilde{J} \tilde{I}})=\text{ det }(- \Gamma _{I J})\text{ det }( V_{ JI})>0$. Finally, using Lemma 24 and Corollary 14, $\tilde{{\mathscr {N}}}$ is robustly non-degenerate.

• 4. Conserved Regulation of a species: We can consider this case as a sequence of two modifications. First, let ${{\hat{{{\mathscr {S}}}}}}= {{\mathscr {S}}}$, and ${{\hat{{{\mathscr {R}}}}}}={{\mathscr {R}}}\cup \{ X_i \leftrightharpoons \emptyset \}$. W.l.o.g, assume that $i=n$. Then, from the previous case it follows that ${{\hat{{{\mathscr {N}}}}}}=({{\hat{{{\mathscr {S}}}}}},{{\hat{{{\mathscr {R}}}}}})$ is robustly non-degenerate. Let ${{\hat{\Gamma }}}$ be the corresponding stoichiometry matrix. Hence, using Lemma 24 and Corollary 14, there exist sets ${{\hat{I}}}, {{\hat{J}}}$ with $\vert {{{\hat{I}}}}\vert =\vert {{{\hat{J}}}}\vert ={\text{ rank }}{{{\hat{\Gamma }}}}$ such that $\text{ det }(-{\hat{\Gamma }}_{{{\hat{I}}} {{\hat{J}}}})\text{ det }({{\hat{V}}}_{{{\hat{J}}} {{\hat{I}}}})>0.$ Next, we define $\tilde{{\mathscr {N}}}$ as follows: $\tilde{{\mathscr {S}}}={\hat{{{\mathscr {S}}}}}\cup \{ X_{n+1}\}$ and $\tilde{{\mathscr {R}}}$ is defined as follows: both $\tilde{{\mathscr {R}}},{{\hat{{{\mathscr {R}}}}}}$ have the same reactions except for $X_n \leftrightharpoons \emptyset$ which is replaced by $X_n \leftrightharpoons X_{n+1}$. Now consider two cases: ${\text{ rank }}(\tilde{\Gamma })={\text{ rank }}({{{\hat{\Gamma }}}})$ and ${\text{ rank }}(\tilde{\Gamma })=1+{\text{ rank }}{({{\hat{\Gamma }}})}$. In first case, using the same argument as in the case of external regulation, $\tilde{{\mathscr {N}}}$ is robustly non-degenerate. We now study the case in which ${\text{ rank }}\tilde{\Gamma }=1+{\text{ rank }}{{{\hat{\Gamma }}}}$. Using a similar argument to the case of external regulation, we can choose ${{\hat{J}}}$ such that $\nu ,\nu +1 \not \in {{\hat{J}}}$, where ${\textbf{R}}_{\nu },{\textbf{R}}_{\nu +1}$ are the reactions $X_{n+1}\rightarrow X_n, X_n \rightarrow X_{n+1}$, respectively. Hence, let $\tilde{I}={{\hat{I}}} \cup \{n+1\}, \tilde{J}={{\hat{J}}}\cup \{\nu \}$. Hence, similar to (25) and (22) we can show that $\text{ det }(-\tilde{\Gamma }_{\tilde{I} \tilde{J}})= \text{ det }(-\tilde{\Gamma }_{{{\hat{I}}} {{\hat{J}}}})$, $\text{ det }({{\hat{V}}}_{ \tilde{J} \tilde{I}})= \text{ det }({{\hat{V}}}_{{{\hat{J}}} {{\hat{I}}} })$. Therefore, $\text{ det }(-\tilde{\Gamma }_{\tilde{I} \tilde{J}})\text{ det }(\tilde{V}_{\tilde{J} \tilde{I}}))>0$. Finally, using Lemma 24 and Corollary 14, $\tilde{{\mathscr {N}}}$ is robustly non-degenerate.

• 5. Adding a catalyst: Let $\tilde{{\mathscr {S}}}={{\mathscr {S}}}\cup \{ X_i^-\}$, and $\tilde{{\mathscr {R}}}$ is defined as in Definition 3, Item 6. It can be seen that this implies that $\tilde{\Gamma }=[\Gamma ^T,-\gamma _i]^T$, where $\gamma _i$ is the ith row of $\Gamma$. Therefore, ${\text{ rank }}\tilde{\Gamma }={\text{ rank }}\Gamma$. Let $\tilde{I}=I, \tilde{J}=J$. Since all the reactions in $\tilde{{\mathscr {N}}}$ are extensions of the corresponding reactions in ${{\mathscr {N}}}$, then the positive term $\text{ det }(-\Gamma _{IJ})\text{ det }(V_{JI})$ is present in the expansion of $\text{ det}_{\tilde{I}}(-\tilde{\Gamma }\tilde{V})$. Therefore, using Lemma 24 we get that $\text{ det}_{\tilde{I}}(-\tilde{\Gamma }\tilde{V})>0$. Hence, using Corollary 14, $\tilde{{\mathscr {N}}}$ is robustly non-degenerate.

• 6. Adding a dimer: Let $\tilde{{\mathscr {S}}}={{\mathscr {S}}}\cup \{ X_i^+\}$, and $\tilde{{\mathscr {R}}}$ is defined as in Definition 3, item 7. It can be seen that this implies that $\tilde{\Gamma }=[\Gamma ^T,\gamma _i]^T$, where $\gamma _i$ is the ith row of $\Gamma$. Therefore, ${\text{ rank }}\tilde{\Gamma }={\text{ rank }}\Gamma$. Using the same argument as in the previous case, we get that $\tilde{{\mathscr {N}}}$ is robustly non-degenerate.

$\square$

1.4 Proof of Theorem 21

Let ${{\mathscr {N}}}=({{\mathscr {S}}},{{\mathscr {R}}})$, and let $\tilde{{\mathscr {N}}}=(\tilde{{\mathscr {S}}},\tilde{{\mathscr {R}}})$ be its elementary modification. For a given reaction ${\textbf{R}}_j$, let ${\mathcal {I}}({\textbf{R}}_j)\subset {{\mathscr {S}}}$ denotes its reactants, while ${\mathcal {O}}({\textbf{R}}_j) \subset {{\mathscr {S}}}$ denotes its products. The statement of the theorem is equivalent to proving that the absence of critical siphons for ${{\mathscr {N}}}$ implies the same for $\tilde{{\mathscr {N}}}$. Pick any $P \subset {{\mathscr {S}}}$. By assumption, P is not a critical siphon for ${{\mathscr {N}}}$. Hence, P is either not a siphon, or it is a trivial siphon. For the first case, using the definition of a siphon, P is not a siphon if and only if the following statement $(\clubsuit )$ holds:

$(\clubsuit )$ $\exists X_k \in P, {\textbf{R}}_k \in {{\mathscr {R}}}$ such that $X_k$ is a product of ${\textbf{R}}_k$ (i.e., $X_k\in {\mathcal {O}}({\textbf{R}}_k)$), and ${\mathcal {I}}({\textbf{R}}_k)\cap P =\emptyset$.

For the second case, if P is a trivial siphon, we assume, w.l.o.g, that it is minimal, i.e., P coincides exactly with the support of a single conservation law. We are ready now to consider the following cases:

• 1. Reversal: By definition, $\tilde{{\mathscr {S}}}={{\mathscr {S}}}$, and for some $j\in \{1,..,\nu \}$ we have $\tilde{{\mathscr {R}}}={{\mathscr {R}}}\cup \{{\textbf{R}}_{-j}\}$. If P is not a siphon for ${{\mathscr {N}}}$, and since $\tilde{{\mathscr {S}}}={{\mathscr {S}}}$ (i.e., no new species added), then the statement $(\clubsuit )$ holds also for $\tilde{{\mathscr {N}}}$. Hence, P is not a siphon for $\tilde{{\mathscr {N}}}$. If P is a trivial siphon for ${{\mathscr {N}}}$, then it is also a trivial siphon for $\tilde{{\mathscr {N}}}$ since addition of a reverse of a reaction does not change the conservation laws of a network. In summary, P is not a critical siphon for $\tilde{{\mathscr {N}}}$. Since $\tilde{{\mathscr {S}}}={{\mathscr {S}}}$, $\tilde{{\mathscr {N}}}$ lacks critical siphons.

• 2. External regulation: $\tilde{{\mathscr {S}}}= {{\mathscr {S}}}$, and $\exists X_i \in {{\mathscr {S}}}$ such that $\tilde{{\mathscr {R}}}={{\mathscr {R}}}\cup \{ X_i \leftrightharpoons \emptyset \}$. Pick any $P \subset {{\mathscr {S}}}$. If P is not a siphon, then the same argument used for the previous modification shows P is not a siphon for $\tilde{{\mathscr {N}}}$. If P is a (minimal) trivial siphon for ${{\mathscr {N}}}$, then either: (A) $X_i \not \in P$ which means that P is a trivial siphon for $\tilde{{\mathscr {N}}}$, or (B) $X_i \in P$ which means that P no longer contains the support of a conservation law for $\tilde{{\mathscr {N}}}$ since $X_i$ has an inflow and is no longer conserved. In summary, P is not a critical siphon for $\tilde{{\mathscr {N}}}$. Since $\tilde{{\mathscr {S}}}={{\mathscr {S}}}$, $\tilde{{\mathscr {N}}}$ lacks critical siphons.

• 3. Conserved regulation: $\tilde{{\mathscr {S}}}= {{\mathscr {S}}}\cup \{X^*\}$, and $\exists X_i \in {{\mathscr {S}}}$ such that $\tilde{{\mathscr {R}}}={{\mathscr {R}}}\cup \{ X_i \leftrightharpoons X^* \}$. If P is not a siphon, then P is not a siphon for the modified network $\tilde{{\mathscr {N}}}$ because the statement $(\clubsuit )$ continues to hold. If we define $\tilde{P}:=P \cup \{X^*\}$, then $\tilde{P}$ is not a siphon since ${\mathcal {I}}({\textbf{R}}_k)\cap \tilde{P} =\emptyset$, i.e., $(\clubsuit )$ holds. If P is a (minimal) trivial siphon for ${{\mathscr {N}}}$, then either: (case A) $X_i \not \in P$ which means that P is a trivial siphon for $\tilde{{\mathscr {N}}}$, or, (case B) $X_i \in P$ which means that P no longer contains the support of a conservation law for $\tilde{{\mathscr {N}}}$ and P is no longer a siphon for $\tilde{{\mathscr {N}}}$. Instead, $P\cup \{X^*\}$ contains the support of a conservation law, and hence, $P\cup \{X^*\}$ is a trivial siphon for $\tilde{{\mathscr {N}}}$. In summary, neither P nor $P\cup \{X^*\}$ are critical siphons for $\tilde{{\mathscr {N}}}$. Since all subsets of $\tilde{{\mathscr {N}}}$ can be represented as P or $P \cup \{X^*\}$ for some $P\subset {{\mathscr {S}}}$, $\tilde{{\mathscr {N}}}$ lacks critical siphons.

• 4. Adding an intermediate: Recall that the network ${{\mathscr {N}}}$ is modified by replacing a reaction of the form ${\textbf{R}}_j=\sum _{i}\alpha _{ij}X_i \rightarrow \sum _i \beta _{ij} X_i$ by two reactions $\tilde{\textbf{R}}_j:= (\sum _{i}\alpha _{ij}X_i \rightarrow X_{n+1})$, and $(\tilde{\textbf{R}}_{\nu +1}:= X_{n+1} \rightarrow \sum _i \beta _{ij} X_i)$. If P is not a siphon for ${{\mathscr {N}}}$, then the statement $(\clubsuit )$ continues to hold for $\tilde{{\mathscr {N}}}$, i.e., P is not a siphon for ${{\mathscr {N}}}$. Next, let $\tilde{P}=P\cup \{X^*\}$. For the sake of contradiction, assume that $\tilde{P}$ is a siphon, this is only possible if $X^*\in {\mathcal {I}}({\textbf{R}}_k)$ (where ${\textbf{R}}_k$ is defined in the statement $(\clubsuit )$). But using our construction, this means that $\mathcal {\textbf{R}}_k={\textbf{R}}_j$. Furthermore, since $\tilde{P}$ is a siphon and it contains $X^*$, one of the reactants of $\tilde{\textbf{R}}_j$ is in $\tilde{P}$. This means that one of the reactants of ${\textbf{R}}_j(={\textbf{R}}_k)$ is in P. But, this contradicts the statement $(\clubsuit )$. Hence, $\tilde{P}$ is not a siphon for $\tilde{{\mathscr {N}}}$. If P is a (minimal) trivial siphon for ${{\mathscr {N}}}$, then it contains the support of a conservation law $d \in {\mathbb {R}}^n_{\ge 0}$. Since P is the support of d, denote $s:=\vert P\vert$, and recall that $\gamma _{ij}$ denotes the (i, j)th entry of $\Gamma$. W.l.o.g, assume that ${{\mathscr {S}}}$ is indexed such that the first s elements coincide with the elements of P. This also implies that $d_1,..,d_s >0$. Hence, $\forall j \in \{1,..,\nu \},~\sum _{i=1}^s d_i \gamma _{i j} =0$. Next, we consider few cases: (Case A) ${\mathcal {I}}({\textbf{R}}_j) \cap P =\emptyset , {\mathcal {O}}({\textbf{R}}_j) \cap P =\emptyset$. We can see that the addition of an intermediate does not change the conservation law; therefore, P is a trivial siphon for $\tilde{{\mathscr {N}}}$, while $P \cup \{X^*\}$ is not a siphon since we assumed that ${\mathcal {I}}({\textbf{R}}_j)\cap P =\emptyset$. (Case B) ${\mathcal {I}}({\textbf{R}}_j) \cap P =\emptyset , {\mathcal {O}}({\textbf{R}}_j) \cap P \ne \emptyset$, or ${\mathcal {I}}({\textbf{R}}_j) \cap P \ne \emptyset , {\mathcal {O}}({\textbf{R}}_j) \cap P =\emptyset$. By AS2, this means either $\sum _{1\le i \le s, \gamma _{ij}>0} d_i \gamma _{i j} =0$ or $\sum _{1\le i \le s, \gamma _{ij}<0} d_i \gamma _{i j} =0$, respectively. Either case contradicts $d_1,..,d_s>0$. (Case C) ${\mathcal {I}}({\textbf{R}}_j) \cap P \ne \emptyset , {\mathcal {O}}({\textbf{R}}_j) \cap P \ne \emptyset .$ Using AS2, we get $\sum _{1\le i \le s, \gamma _{ij}<0} d_i \gamma _{i j}=\sum _{1\le i \le s, \gamma _{ij}>0} d_i \gamma _{i j}:=\xi >0$. Therefore, $\tilde{P}=P\cup \{X^*\}$ is a trivial siphon for $\tilde{{\mathscr {N}}}$ with the conservation law ${{\tilde{d}}}=[d^T \xi ]^T$, while P is not a siphon for $\tilde{{\mathscr {N}}}$. In summary, neither P nor $P\cup \{X^*\}$ are critical siphons for $\tilde{{\mathscr {N}}}$. Since all subsets of $\tilde{{\mathscr {N}}}$ can be represented as P or $P \cup \{X^*\}$ for some $P\subset {{\mathscr {S}}}$, $\tilde{{\mathscr {N}}}$ lacks critical siphons.

• 5. Adding a dimer: Let $\tilde{{\mathscr {S}}}={{\mathscr {S}}}\cup \{ X_i^+\}$, and $\tilde{{\mathscr {R}}}$ is defined as in Definition 3, item 7. If P is not a siphon, this means that the statement $(\clubsuit )$ holds. Since $X_i^+$ shares the same input and output reactions with $X_i$, the statement $(\clubsuit )$ continues to hold. Hence, P is not a siphon for $\tilde{{\mathscr {N}}}$. If P is a trivial siphon for ${{\mathscr {N}}}$ then it will be a trivial siphon for $\tilde{{\mathscr {N}}}$ since adding a dimer preserves the existing conservation laws of ${{\mathscr {N}}}$. Now, let us consider a set of the form $\tilde{P}:=P\cup \{X_i^*\},P\subset {{\mathscr {S}}}$. For the sake of contradiction, assume that $\tilde{P}$ is a critical siphon for $\tilde{{\mathscr {N}}}$. Let ${{\hat{P}}}:= (\tilde{P}/\{X_i^+\})\cup \{X_i\}$. Since $X_i$ have the same reactants and products as $X_i^+$ by construction, then ${{\hat{P}}}$ is a siphon for ${{\mathscr {N}}}$ and it does not contain the support of a conservation law. Hence, ${{\hat{P}}}$ is a critical siphon for ${{\mathscr {N}}}$ which contradicts our assumption. Since all subsets of $\tilde{{\mathscr {N}}}$ can be represented as P or $P \cup \{X^*\}$ for some $P\subset {{\mathscr {S}}}$, $\tilde{{\mathscr {N}}}$ lacks critical siphons.

1.5 Proof of Theorem 22-1

Part 1 of Theorem 22 follows from Theorem 21 except for the case of adding a catalyst which is proved next. Assume that $\tilde{{\mathscr {N}}}$ is a modification of a linear network ${{\mathscr {N}}}$ by adding a catalyst as described in Definition 3. Let $P \subset {{\mathscr {S}}}$. So P is either not a siphon or a trivial siphon. Assume first that P is not a siphon, hence the statement $(\clubsuit )$ in proof of Theorem 22 holds. Since the products and reactants of reactions in $\tilde{{\mathscr {N}}}$ contain their counterparts in ${{\mathscr {N}}}$ then P is not siphon for $\tilde{{\mathscr {N}}}$. Consider $\tilde{P}=P\cup \{X_i^-\}$. We consider two cases: (case A) $X_i \not \in \tilde{P}$. Note both ${{\mathscr {N}}}$ and $\tilde{{\mathscr {N}}}$ satisfy AS1. Therefore, $\tilde{{\mathscr {R}}}$ must contain at least one reaction $\tilde{\textbf{R}}_j$ with $X_i^- \in {\mathcal {O}}(\tilde{\textbf{R}}_j)$. In order for $\tilde{P}$ to be a siphon, there must exist $X^* \in {\mathcal {I}}(\tilde{\textbf{R}}_j) \cap \tilde{P}$. But since ${{\mathscr {N}}}$ is linear, ${\textbf{R}}_j$ has at most one reactant and it must be $X_i$ (since by construction, input reactions of $X_i^-$ are output reactions of $X_i$). This also implies that $\tilde{\textbf{R}}_j$ has at most one reactant, and hence $X_i=X^*\in \tilde{P}$ which is a contradiction. (case B) $X_i \in \tilde{P}$; hence, $\tilde{P}$ contains the support of the conservation law $X_i + X_i^-=\text{ constant }.$ Therefore, $\tilde{P}$ is a trivial siphon, and it is not critical.

Second, let us assume that P is a trivial siphon for ${{\mathscr {N}}}$, then P and $P\cup \{X_i^-\}$ are trivial siphons for $\tilde{{\mathscr {N}}}$ since adding a catalyst preserves the existing conservation laws of ${{\mathscr {N}}}$. Since all subsets of $\tilde{{\mathscr {N}}}$ can be represented as P or $P \cup \{X^*\}$ for some $P\subset {{\mathscr {S}}}$, $\tilde{{\mathscr {N}}}$ lacks critical siphons.