Many biological systems are known for the ability to operate precisely and consistently subject to potentially large disruptions and uncertainties [1][2][3][4][5]. Examples are homeostasis, understood as the maintenance of a desired steady state (perhaps associated to an observable phenotype) against the variability of in-vivo concentrations of biochemical species, or a consistent dynamical behavior in the face of environmental variations which change the speed of reactions.
The vaguely defined term "robustness" is often used to refer to this consistency of behavior under perturbations. The present work deals with such notions of "biological robustness", as well with a "robustness of analysis" notion in which conclusions can be drawn despite inaccurate mathematical models. Models of core processes in cells are typically biochemical reaction networks. This includes binding and unbinding, production and decay of proteins, regulation of transcription and translation, metabolic pathways, and signal transduction [6]. However, in contrast to engineered chemical systems, biology poses particular challenges. On the one hand, the reactants and the products in such interactions are frequently known, and hence the species-reaction graph is available. On the other hand, the exact form and parameters (i.e., kinetics) that determine the speed of transformation of reactants into products are often unknown. This lack of information is a barrier to the construction of complete mathematical models of biochemical dynamics. Even if the kinetics are exactly known at a specific point in time, they are influenced by environmental factors and hence they can change. Hence, the ability to draw conclusions regarding the qualitative behavior of such networks without knowledge of their kinetics is highly relevant, and has been advocated under the banner of "complex biology without parameters" [4]. But is such a goal realistic? It is known that the long-term qualitative behavior of a nonlinear system can be critically dependent on parameters, a phenomenon known as bifurcation. This fundamental difficulty led to statements such as Glass and Kauffman's 1973 assertion that "it has proved impossible to develop general techniques which may be applied to find the asymptotic behavior of complex chemical systems" [7].
Notwithstanding such difficulties, many classes of reaction networks are observed to have a "well-behaved" qualitative long-term dynamical behavior for wide ranges of parameters and various types of nonlinearities. This means specifically in our context that such networks do not have the potential for exhibiting complex steady-state phenotypes such as multiple-steady states (e.g., toggle switches), oscillations (e.g., repressilator), or chaos. Their typical behavior is that the concentrations eventually settle into a unique steady state (called an attractor) for any initial condition (with fixed total substrate, gene and enzyme concentrations). Hence, we call them structurally attractive. The relevant biological phenotype for such networks is the unique attractor, which is mathematically represented by the concentrations of the biochemical species at steady state. Discerning such networks is not generally trivial. For instance, within the class of post-translational modification (PTM) cycles, some cascades are "structurally attractive" but others can exhibit oscillations and multistability [8]. Fig 1 illustrates the typical behavior of an attractive network vs a multistable network for two PTM cycles that have been proposed as models for double phosphorylation. We will study PTM cycles in detail later in the paper.
In the terminology of dynamical and control system theories, the defining feature of an attractive network is that it can only exhibit global point attractors (i.e., unique globally asymptotically stable steady states). The classical way to certify stability is by exhibiting an appropriate energy-like function, commonly referred to as a Lyapunov function [9,10]. Existence of such a function provides many guarantees on qualitative behavior, including notably the fact that its sub-level sets act as trapping sets for trajectories [11]. Furthermore, they allow the development of a systematic study of model uncertainties and response to disturbances [9,10]. However, it is notoriously difficult to find such functions for nonlinear systems due to the lack of general constructive techniques.
The search of Lyapunov functions for nonlinear reaction networks can be traced back to Boltzmann's H-Theorem [12], which applies only to the restrictive subclass of detailed-balanced networks. Wei [13] in 1962 postulated that all chemical systems should satisfy an "axiom of convergence" and there shall exist a suitable Lyapunov function. Perhaps the most striking success in this line of thought was the development of the Horn-Jackson-Feinberg (HJF) theory of complex-balanced networks [14][15][16][17] in the early 1970s, which relies on using the sum of all the chemical pseudo-energies stored in species as a Lyapunov function. When Distinct qualitative behaviors for two models of a double PTM. This is illustrated by the time series plots for the double phosphorylated substrate with randomized initial conditions for fixed total substrate and enzyme concentrations. (a) the processive mechanism exhibits a unique global attractor, (b) a distributive mechanism exhibits multistability for some parameters. See networks (11), (12) and the accompanying discussion. The parameters are given in S1 Text- §6. https://doi.org/10.1371/journal.pcbi.1007681.g001 specific graphical conditions are satisfied, complex-balancing is guaranteed for all kinetic constants. Global stability can be proven in certain cases [18,19]. Despite the elegance and theoretical appeal of the method, the assumptions needed for its applicability are restrictive, and are not widely satisfied in biological models. For example, many basic motifs (e.g., transcription/translation and enzymatic reactions) are not complex balanced. Furthermore, HJF theory assumes, although with some exceptions, that the reaction kinetics are Mass-Action. It has been argued that this assumption "is not based on fundamental laws" and is merely "good phenomenology" [20]. These laws are usually justified by the intuitive image of colliding molecules. However, this is often not the right level of analysis for biological modeling, where alternative kinetics such as Michaelis-Menten and Hill kinetics are used in situations involving multiple time scales [21].
Beside complex-balanced networks, a few additional classes of attractive networks have been identified. These include mono-molecular networks, which can be handled within the framework of compartmental systems using a Lyapunov function [22,23]. More recently, global convergence has been shown for another class of networks via the concept of monotonicity without supplying a Lyapunov function [24] where sufficient graphical conditions have been given.
In previous work [25][26][27], two of the authors proposed a direct approach to the problem, introducing the class of piecewise linear-in-rates functions, which act as Lyapunov functions regardless of the specific form of the reaction nonlinearities or kinetic constants. They guarantee the uniqueness of steady states and global stability under mild additional conditions.
In this work, the results from [25][26][27] are generalized in several directions, theoretically, computationally, and in terms of biological applications. First, we propose a general characterization of "structurally attractive" networks. We require the existence of a universal ratedependent function, which we call a Robust Lyapunov Function (RLF), that is a Lyapunov function for any choice of the kinetics. We proceed to propose a computational framework for finding such functions. To this end, the dynamics of the network are embedded in a linear convex cone. The extremals of this cone are a set of rank-one matrices that derive from the stoichiometry of the network. If a common Lyapunov function exists for the extremals, then it can be used to construct an RLF and the network is certified to be attractive. In the special case that kinetics are mono-or bimolecular, the RLF is piecewise linear or piecewise quadratic on species, respectively.
Computationally, we complement previous reaction network toolboxes [28,29] and we provide a Lyapunov-Enabled Analysis of Reaction Networks (LEARN) toolbox to implement the results on any given network by searching for an RLF and checking the appropriate conditions via four main methods: a graphical algorithm, a linear program, an iterative procedure, and a semi-definite program. Additionally, LEARN checks for conditions that rule out the existence of an RLF.
We then proceed to carry out an extensive discussion of biochemical networks to show the applicability of our results. Foundational studies in systems biology [6] have revealed that biochemical networks have many common "motifs". We show that our results form a basis for the understanding of the behavior of a large class of networks of various degrees of complexity. They may be applied to study basic motifs such as binding/unbinding, three-body binding, transcription and translation networks, and enzymatic reactions. Most cellular signalling involves PTMs as building blocks, and their malfunction is frequent in diseases such as cancer and Alzheimer [30,31]. Hence, we study in detail PTMs cascades, ERK signalling, and phosphotransfer and phosphorelay networks. In addition, we study important biological networks such as T-cell kinetic proofreading, and the Ribosome Flow Model. We show that our Lyapunov functions can be used to construct safety sets and perform dynamic flux analysis. Many of the networks studied are not amenable to the previously-mentioned analysis techniques, HJF theory in particular. A comparison with other methods in included in the Discussion (see Table 1). In particular, our results include the class of monomolecular networks treated in [22,23], and it applies all the biochemical networks studied in [32], [24], [33]. A preliminary version of a subset of these results were presented in conferences [34], [35].
Theoretically, our results connect with a corpus of previous literature. We show that the RLFs can be formulated in different coordinates, and how this relates to the ones proposed in [34], [36]. Also, the approach makes contact with the notions of structural injectivity [37][38][39][40], structural persistence [41], and uncertain linear systems [42][43][44].
The paper has been written for a diverse readership, and has been structured accordingly. Readers who are interested in the general concepts, the biological applications, and the software package only need to consult the Introduction, the Results, and LEARN's accompanying manual (SI §7). Users can apply the results by supplying the list of reactions encoded as a stoichiometry matrix as an input to LEARN's main subroutine for a report of results. Readers who are also interested in the technical mathematical details can consult the Methods section.
Since LEARN guarantees that a certain mechanism cannot admit multistability, oscillation, or chaos, it can be used to distinguish competing biochemical reaction networks at the modeling step. We give an example of this when discussing processive vs distributive post-translational cycles.
LEARN can be compared to other results in the literature as shown in Table 1.
A list of reactions can be abstracted mathematically into the framework of Chemical Reaction Networks (CRNs). A CRN consists of a set of species S ¼ fX 1 ; ::; X n g and a set of reactions R ¼ fR [50]). The graph corresponding to the Comparison with other methods in the literature. The row that corresponds to "admissible kinetics" asks about the functional form of the reaction rates for which the method is applicable. "Global attractor" asks whether the method is able to provide guarantees for the global convergence to an attractor. "Uniqueness with i/o perturbations" asks whether the method can guarantee uniqueness of steady states with respect to arbitrary addition of inflows and outflows to the network (i.e., "homogeneous CFSTR" in the terminology of [45]). Rows that correspond to "PTM cycle" and "Kinetic proofreading" ask whether the method can tackle the networks ( 9) and ( 15), respectively. We have picked these two networks as non-trivial examples that are pertinent to systems biology. The question of a global attractor for HJF-type networks is marked by an asterisk ( � ) since a proof has been proposed in a preprint [46] but is not formally published yet. PTM cycle is given in Fig 2c). The stoichiometry matrix Γ becomes the incidence matrix of the graph [51].
As we are interested in studying the long-term dynamical behavior, a concentration x i � 0, i = 1, .., n is assigned to each species. Hence, the concentration vector at time t is x(t) = [x 1 (t), . . ., x n (t)] T . A reaction rate (or flux) R j (x), j = 1, .., ν is assigned to each reaction. The reaction rate vector is R(x) = [R 1 (x), . . ., R ν (x)] T . The time-evolution of the concentration vector is given by the standard ordinary differential equation (ODE) given as [52]:
Biochemical networks usually contain conserved quantities (i.e., moieties) such as the total amount of enzymes, substrates, ribosomes, RNA polymerase, etc. For each conserved quantity, there exists a nonnegative vector d such that d T Γ = 0, and d is called a conservation law. If every species is supported in at least one conservation law the network is said to be conservative. For example, the PTM cycle in Fig 2 is conservative with three conservation laws c 1 + c 2 + x + y = [X] total , e + c 1 = [E] total , and f + c 2 = [F] total , which are the total amounts of the substrate and the two enzymes, respectively, and they stay constant throughout the reaction. Hence, claims of global stability and uniqueness of steady states are relative to the conserved quantities. A set of concentrations that shares the same conserved quantities is called a stoichiometric class.
For the PTM cycle, the ODE is given in Fig 2b). We do not assume that the reaction rates have a specific functional form such as Mass-Action. We only assume that the rates are monotone, meaning that as the concentration of reactants increases, the rate of the reaction increases (see Methods). This can be interpreted as enforcing a specific sign pattern on the partial derivatives of R. This means that all the entries of the Jacobian matrix of R (i.e., @R/@x), are either zero or non-negative. For the PTM cycle, Fig 2d) illustrates our assumptions on the reaction rates encoded in terms of the Jacobian matrix. Such reactions include all common reaction rates such as Mass-Action, Michaelis-Menten, Hill, etc.
Despite its application relevance, establishing the long-term behavior of the PTM cycle in Fig 2 was an open problem till the 2000s. HJF's theory cannot be used for deciding stability since the PTM cycle is a non-zero deficiency network. In 2008, this problem was tackled via monotonicity techniques [24,32], but no Lyapunov function has been provided. As a motivation, we study the same cycle using our proposed method. An intuitive way to approach its analysis is to consider the central loop in Fig 2, and then study the sum of absolute rate differences along it. This can be loosely motivated by considering the reactions rates as potentials and the concentration of species as charges, and noting that the difference of "potentials" causes the concentration of species to change via the flow of a "current". Hence, we define the ith current as the rate of change of the concentration of the ith species. Thus, we consider the weighted sum of currents P i w i j _ x i j as a candidate Lyapunov function. It can also be written as follows:
which is a piecewise linear-in-rates function. In order to verify whether this is indeed a Lyapunov function, we can analyze it region-wise to check that it decreases along trajectories. Consider for instance the region W ¼ fR 1 ðxÞ � R 2 ðxÞ � R 3 ðxÞ � R 4 ðxÞg. The candidate V simplifies to the difference of "potentials" across the substrate S:
To evaluate _ V , we need the signs of the "currents" _ s; _ e; _ c 2 . In our example, we can use the inequalities defining W so that the signs can be read from the graph as follows: _ s; _ e < 0 and _ c 2 > 0. By noting that these signs are matched to the coefficients of R(x) in (3), and since @R/ @x is nonnegative, we can write the following inequality in W:
where the sign of the rate of change of each concentration is indicated above it. Therefore, sgn _ V can be determined conclusively without knowing the kinetics. In fact, this can be repeated for all regions to conclude that V is non-increasing along all possible trajectories of (1). (See the Results section for further analysis).
The lesson that can be drawn from this example is that a robust analysis of reaction networks can be carried out by considering candidate Lyapunov functions of the form Ṽ ðRðxÞÞ that vanish exactly on the steady state set, i.e. the set {x|ΓR(x) = 0}. This approach does not require the computation of the actual steady state.
The motivating example has shown that we can have a Lyapunov function Ṽ ðRðxÞÞ that decreases along trajectories for any monotone kinetics R. Hence for a given network ðS ; RÞ we will be looking for a function Ṽ : R n ! R �0 that vanishes only on the set of steady states, i.e Ṽ ðrÞ ¼ 0 if and only if rG ¼ 0: Furthermore, VðxÞ ¼ Ṽ ðRðxÞÞ needs to be nonincreasing along the trajectories of (1), i.e it must satisfy: _ Ṽ ðRðxÞÞ ≔ ð@ Ṽ =@RÞð@R=@xÞGRðxÞ � 0; for all x and for all R admissible: ð4Þ
If such a function exists then we call it a Robust Lyapunov Function (RLF), and the network is called structurally attractive. Mathematically, the RLF needs only to be locally Lipschitz and the derivative is defined in the sense of Dini's (see Methods).
Example (cont'd). For the PTM cycle (Fig 2)
The above definition of an RLF does not offer a constructive way for finding one or for checking a candidate. Our first result is to give a characterization of RLF in terms of a set of rankone linear systems, each of which corresponds to a reaction-reactant pair. The set of all such pairs is P ≔ fðj; iÞjX i participates in the reaction R j g. Let s be total number of such pairs. Then,
., s where {γ 1 , .., γ n } are the rows of Γ and {e 1 , .., e ν } are columns of the ν × ν identity matrix.
The matrices Q 1 , .., Q s will serve as system matrices for s linear systems and also as extremals of a linear convex cone. We show (see Methods) that (@R/@x)Γ 2 cone(Q 1 , .., Q s ) = {∑ ℓ ρ ℓ Q ℓ |ρ ℓ � 0}. We will be looking for a function Ṽ that acts as a common Lyapunov function for these linear systems and satisfies frj Ṽ ðrÞ ¼ 0g ¼ \ s '¼1 ker Q ' (see Methods). We are ready to state the main result of this section. (See Methods). Theorem 1. Given ðS ; RÞ. Let (1) be the associated ODE. A function
is an RLF for the reaction network ðS ; RÞ.
The characterization provided in Theorem 1 can be used for devising computational algorithms that search for an RLF. In Methods, we present several algorithms for constructing piecewise linear (PWL) or piecewise quadratic RLFs. In order to simplify the presentation, we will be only looking for convex piecewise linear RLFs in our study of biochemical networks. This means looking for vectors c 1 ; :::; c m 2 R n (for some positive integer m) such that Ṽ is an RLF where: [22] is a special case with ξ = 1. The vector ξ can be found by linear programming using a special case of Theorem 2 (see Methods). Note that the function (2) discussed in the motivating example has the form ( 6) above.
Max-Min RLFs. These are functions of the form:
where R consists of reaction rates or the difference between forward and backward rates of a reaction. Unlike SoC RLFs which keep track of the reaction rate differences across each species, the Max-Min RLF keeps track of the maximal reaction rate difference across the whole network at each time. We provide a full graphical characterization of the class of networks that admit Max-Min RLFs (which we call M-networks). (see Methods, Theorem 4). Alternative forms. In Methods, we give conditions on a function V such that V ðx x e Þ (where x e is a steady state) is a Lyapunov function for any admissible R. We call V a concentration-dependent RLF. We show that Ṽ ðrÞ Theorem 11). These PWL functions relate to the ones proposed in [34,36]. Note, however, that V ðx x e Þ is a Lyapunov function only in the stoichiometric class that contains x e .
Properties of RLFs. In [27], some properties of networks admitting PWL RLFs have been established and they can serve as necessary condition tests. In Methods, we provide two additional properties, namely testing robust non-degeneracy and the absence of critical siphons. These conditions are implemented in LEARN.
The existence of an RLF implies that the qualitative long-term behavior of a network is highly constrained. Hence, an important issue is whether this theory is sufficiently relevant to biomolecular applications. We will show in the remainder of the Results section that this class of networks constitutes a rich and relevant class. It includes basic motifs, modules, and larger networks and cascades in molecular biology. For most of these networks, the HJF Lyapunov function [14] does not apply. And if it applies, it is only valid with Mass-Action kinetics (or a generalization [18]) and it does not confer the same powerful conclusions offered by our theory. Many of the networks discussed in the remainder of this paper are qualitatively analyzed for the first time and most of them had no Lyapunov functions known for them. For all the subsequent networks the following statement holds: if a positive steady state exists, then it is unique and globally asymptotically stable relative to its stoichiometric class.
In this subsection, several biochemical networks are presented. They are fairly simple and all of them can be analyzed using HJF theory in the case of Mass-Action kinetics. However, they are presented here to show that the properties that our theory requires are obeyed by the basic biochemical motifs, which establishes its applicability and generality. Furthermore, we offer an intuitive window to the meaning of RLFs and how our graphical conditions apply.
which can represent an enzyme binding to a substrate. The corresponding RLF can be found easily using Theorem 4 and is given by:
Both the Max-Min and the SoC RLFs coincide in this case.
By considering the irreversible subnetwork 0 ! E, 0 ! X, X + E ! XE, XE ! 0, a Max-Min RLF can be found using Theorem 4 and is given by (7) where
Cooperative binding reaction. The following reactions (depicted in Fig 3c ) represent the situation where n enzyme molecules E need to bind to each other to react to X:
The case n = 2 is called dimerization. The corresponding RLF can be found using Theorem 4 and R is given by (8). The irreversible subnetwork for which Theorem 4 was applied is 0
The corresponding RLF can be found using Theorem 4 and R is given by ( 8). The irreversible subnetwork for which Theorem 4 was applied is 0
We have applied our techniques to the dynamics of simple binding which can be analyzed easily using various known ways. However, it is often the case that two compounds X, Y cannot bind unless a bridging molecule E allows them to bind, forming a ternary complex. This is known as three-body binding [53] and it is ubiquitous in biology. Examples include T-cell receptors interaction with bacterial toxins [54], coagulation [55], and multi-enzyme supramolecular assembly [56]. The same reaction network also models the binding of two different transcription factors into a promoter with a double binding site. Despite its simplicity, the steady-state analysis of the equilibria has been subject of great interest [53]. Stability cannot be decided via HJF theory, and it has not been studied before to our knowledge.
The network can be depicted in Fig 4, and is given by eight reactions as follows:
The network is an M-network and the corresponding irreversible subnetwork has the reactions {R 1 , R -2 , R -3 , R 4 }. Hence we apply Theorem 4 to have an RLF of the form (7) where
It can be concluded that there exists a unique steady state in each stoichiometric class and it is globally asymptotically stable.
Transcription and translation are the first two essential steps in the central dogma of molecular biology, and hence they are of utmost importance in the analysis of gene regulatory networks.
Transcription network. Fig 5a ) shows the transcription network which describes the production of mRNA from DNA using the RNA polymerase [57]:
This model explicitly accounts for the concentration of RNA polymerase and hence it extends to situations in which RNA polymerase is not abundant.
Applying Theorem 4, the RLF (7) can be used with R ¼ fR 1 R 1 ; R 2 ; R 3 g. Alternatively, Theorem 2 can be used, and the Lyapunov function found can be written as:
, where the species are ordered as RNAP, DNA, RD, mRNA. Note this network has deficiency one, hence no information regarding stability can be inferred from HJF theory. Furthermore, the procedure proposed in [36] has been reported not to work for the network above.
Translation network with a leak. Fig 5b ) shows the translation network which describes the production of a protein from mRNA via ribosomes [57]. The leaking of the Ribosome-mRNA complex into the pool of ribosomes is also modeled. In order to make the model more general, we also explicitly account for the concentrations of ribosomes. This is relevant to situations in which ribosomes are not highly abundant which can occur naturally [58,59] or in synthetic circuits [60]. The network can be written as
Note that the flux corresponding to reaction R 4 vanishes at steady state which implies that the species mRNA:Ribo vanishes at any steady state. Note also that the dynamics of other species are independent of the dynamics of P. Hence, the network can be considered as a
Note that Ṽ is neither SoC nor Max-Min. The second network can be analyzed using this Lyapunov function:
Overall stability is established for the cascade using standard techniques [61].
Basic activation motif. Fig 6a ) represents the basic enzymatic reaction where an enzyme E binds to a substrate S to produce S + as follows [48]:
Theorem 3 can be used. The resulting Lyapunov function is:
Although this network has deficiency zero, it is not weakly reversible. This implies that the steady states belong to the boundary, and HJF theory does not offer any information regarding stability.
Enzymatic activation cycle. In order to close the cycle of the activation motif, Fig 6c) depicts the activation of a protein P by an enzyme E, and then the activated protein decays back to its inactive state. The list of reactions is given as [62]:
Theorem 2 gives the following SoC RLF:
and both Theorems 3 and 4 give RLFs also. This network has deficiency one; the deficiency one algorithm [17] excludes the existence of multiple steady states with Mass-Action kinetics. No information regarding stability can be inferred in that context from HJF theory. Furthermore, the decay reaction R 3 usually models fast dephosphorylation which has a Michaelis-Menten kinetics, which is not allowed in [17].
The full PTM cycle. A simplified version of the enzymatic futile cycle has already been used as a motivating example in Fig 2 . It differs from the preceding network by explicitly modeling the dephosphorylation step. The following describes the complete model [48,49]:
For instance, S represents the base substrate, E is called a kinase which adds a phosphate group to S to produce S + . This process is called phosphorylation. The dephosphorylation reaction is achieved by a phosphatase F that removes the phosphate group from S + to produce S.
Theorem 4 can be used to find the RLF (7) where
Alternatively, Theorems 3 yields the SoC RLF:
Both SoC and Max-Min RLFs have an intuitive meaning in terms of the reaction graphs of the networks. The first is the difference between the fastest and the slowest reactions, and the second is the sum of currents (rates of change of concentrations). Since the deficiency of the network is one, stability cannot be inferred from HJF theory.
Energy-constrained PTM cycle. Basic Motif. Madhani [63] presents this biochemical example of adding a phosphate group to a protein using a kinase. ATP is not assumed to be abundant and its dynamics are explicitly modeled. The reaction network is depicted in black in Fig 7a), which can be written as:
where K is the kinase, ATP is the adenosine triphosphate, ADP is the Adenosine diphosphate, and P + is the phosphorylated protein. Reactions R 3 , R 4 are not supported in the kernel of the stoichiometry matrix, which implies that the species PAK, P + A -K vanish at any steady state point.
Applying Theorem 3, one can get the following RLF function:
VðxÞ ¼ max fjR 1 ðxÞ R 1 ðxÞj; jR 2 ðxÞ R 2 ðxÞj; R 3 ðxÞ; R 4 ðxÞ; jR 5 ðxÞ R 5 ðxÞjg:
Energy constrained PTM cycle. In order to have a full cycle, the model can include the following two reactions:
, where ADP is converted to ATP by other cellular processes and is modeled as a single step, and P + decays to its original state P spontaneously or chemically [64]. The reaction network is depicted in Fig 7a).
The full network is an M network, and it has the RLF (7) with
The network is not complex-balanced and HJF theory is not applicable. The dynamics of this network have not been analyzed before per our knowledge.
Full energy-constrained PTM cycle. The dephosphorylation step can be modeled fully and is depicted in Fig 7b). This is the energy-constrained analog of Fig 6c). The network is also an M-network and it admits an RLF of the form (7). The list of reactions have not been included for the sake of brevity.
The post-translation modification (PTM) cycle (e.g, phosphorylation-dephosphorylation cycle [48,49]) has been analyzed in the previous section. This kind of cycle appears frequently in biochemical networks, and can be interconnected in several ways; we discuss some here. For recent reviews see [65,66].
A multisite PTM with distinct enzymes. It is known that a single protein can have up to different 100 different PTM sites [65] and it can undergo different PTM cycles such as phosphorylation, acetylation and methylation [67,68]. Each of these cycles has its own enzymes.
Hence, we consider a cascade of n PTM cycles as shown in Fig 8a ) where n is any integer greater than zero. For instance, the associated reaction network for the case n = 2 is given as:
The network is not an M-network and hence Theorem 4 is not applicable. However, using Theorem 2 it can be shown that a SoC RLF for the n cascade exists and can be represented as VðxÞ ¼k diagðxÞ _ xk 1 with ξ = [2, 2, . . .., 2, 1, 1, . . ., 1] T with the ordering given as X 0 , . . ., X n , E 0 , E 1 , . . ., F n-1 X n . HJF theory will not apply since this network has deficiency n. Also, monotonicity-based results [24] do not apply, since the network is not cooperative in reaction coordinates. In fact, the long-term behavior of this cascade has not been studied before to our knowledge. It follows that for any n the network has a unique globally asymptotic stable steady state in any stoichiometric class (i.e., with respect to fixed total amounts for the enzymes and the substrate).
Multiple PTM cycle with a processive mechanism. Proteins can undergo different PTMs, but they also can undergo a multisite PTM. For instance, a phosphate group can be added to multiple sites on the protein [69]. Multisite phosphorylation can be processive [70] or distributive [71]. The reaction network can be written as [33]
It can be noticed that for every n, the network satisfies the graphical conditions of Theorem 4. Therefore, an RLF is (7) where
Energy-constrained processive cycle. The ATP and ADP expenditure can be accounted for in the processive cycle similar to the model presented in Fig 7b). The new network will remain an M-network and Theorem 4 can be applied. Details are omitted for brevity.
A generalized processive cycle. An "all-encompassing" processive cycle has been studied in [8] which allows multiple enzymes and is depicted in Fig 8c . It takes the following form:
. . .
This network is also an M network and it satisfies the results of Theorem 4. Hence, the Lyapunov function ( 7) can be used.
Both networks above have been studied in [8,33] by establishing monotonicity in reaction coordinates. Such techniques require checking persistence a priori and do not provide Lyapunov functions. Furthermore, our results have the advantage of providing an "all-encompassing" general framework that includes many of these individually studied networks in addition to new ones.
Distinguishing between processive and distributive mechanisms. Fig 8d ) depicts a double futile cycle with a distributive mechanism [71,72], which is described by the following set of reactions:
It can be verified that the network violates the P 0 necessary condition (for the minor corresponding to X 0 , X 1 , X 2 , E, FX 1 , EX 1 ). Hence, a PWL RLF does not exist [27]. Indeed, the above network is known to admit multi-stability for some parameter choices as shown in Fig 1.
Hence, our results can be used to compare between distributive and processive mechanisms as viable models for the first stage in the MAPK cascade. Since the latter has been observed experimentally to accommodate multiple non-degenerate steady states, the processive mechanism cannot be a model.
Phosphotransfer is a covalent modification in which a histidine kinase gives the phosphate group to a response regulator and it is the core motif in a two-component signaling systems [75]. Phosphotransfer cascades are called phosphorelays [76,77]. Phosphotransfer motif. An example is the envZ/ompR signaling system for regulating osmolarity in bacteria such as E. Coli [78]. The core motif can be described by the following set of reactions [79]:
where the "+" superscript refers to a phosphorylated substrate. For instance, Z + is the phosphorylated EnvZ protein, while X is the ompR protein.
The proteins Z, X + can also be phosphorylated and dephosphorylated by other reactions. Fig 9a) presents a network where those other reactions are modeled as a single step:
where R 3 (which phosphorylates Z) can be monotonically dependent on external signals such as osmolarity in the envZ/OmpR network.
It can be noticed that Theorem 4 is applicable and ( 7) is an RLF with
Phosphotransfer with enzymes. A more elaborate model can into account the phosphorylation/dephosphorylation of proteins Z, X + in terms of other enzymes. Hence, reactions (13) can be replaced by the following:
as depicted in Fig 9b . Similarly, ( 7) is an RLF with
A phosphorelay is a cascade of several phosphotransfers. It appears ubiquitously in many organisms. For example, the KinA-Spo0F-Spo0B-Spo0A cascade in Bacillus subtilis [80] and the Sln1p-Ypd1p-Ssk1p cascade in yeast [81].
Fig 9c depicts the cascade which is given by:
where the first kinase is phosphorylated by some constant external signal, and X þ n is the response regulator.
The network is still an M-network and conditions of Theorem 4 apply by mere inspection of the graph. Hence a function of the form ( 7) is a Lyapunov function. Enzymatic activation/ deactivation of X 1 ; X þ n , respectively, can also be added (analogously to Fig 9b ) and the result will continue to hold. Note that the same applies to the more general model presented in [82] also. We omitted the details for brevity.
Note that none of the phosphotransfer networks is complex-balanced and hence HJF theory is not applicable.
In 1974, Hopfield [83] proposed the kinetic proofreading model in protein synthesis and DNA replication. Subsequently, McKeithan [84] proposed a network containing a ligand, which is a peptide-major histocompatibility complex M, binding to a T-cell receptor; the receptor-ligand complex undergoes several reactions to reach the final complex C N . The chain of reactions enhances the recognition and hence it is called a kinetic proofreading process. Fig 10a) depicts the reaction network, which is given by the following set of reactions:
Applying Theorem 2, it can be shown that for any N � 1, the network admits a SoC RLF of the form V N ðxÞ ¼ k diagð½1; 1; 2; 2; ::; 2� T Þ _ x k 1 , where the species are ordered as T, L, C 0 , C 1 , . . ., C N . Note that this network does not meet the graphical requirements of Theorem 4 since it is not an M network. The monotone-systems approach proposed in [24] is not applicable here since the system is not cooperative in reaction coordinates.
Nevertheless, this is one of the few networks, considered so far, which is complex-balanced. The work [18] showed that this network is weakly reversible and that it has zero-deficiency; therefore any positive steady state is unique relative to the interior and is locally asymptotically stable. In order to infer global stability, it was necessary to compute the steady states explicitly to preclude a boundary steady state stoichiometrically compatible with a positive steady state. In comparison, our approach is more powerful, since the former approach is limited to generalized Mass-Action kinetics, and cannot infer global stability directly.
Fig 10b depicts the network describing the effect of the so called Raf Kinase Inhibitor Protein (RKIP) on the Extracellular Regulated Kinase (ERK) signaling pathway as per the model given in [85]. It can be described using the network:
where K is the RKIP, E is the ERK Kinase, P is the RKIP phosphatase, and M is the phosphorylated MAPK/ERK Kinase, and the plus superscript means that the molecule is phosphorylated.
The network is an M-network and the requirements of Theorem 4 are satisfied. Hence, ( 7) is an RLF with R ¼ fR k R k ; k ¼ 1; ::; ng, and R -k (x) :� 0 if R k is irreversible. Note that this network is of deficiency one, hence stability cannot be inferred by HJF theory. Nevertheless, monotonicity-based analysis can be applied [24] which utilizes cooperativity in reaction coordinates. Refer to the Discussion for a detailed comparison to monotonicity techniques.
Finally, we show that our techniques' applications in molecular biology are not limited to classical biochemical networks. A translation elongation process involves ribosomes travelling down an mRNA, readings codons and translating amino-acid chains via recruited tRNAs. A conventional stochastic model is the Totally Asymmetric Simple Exclusion Process [86]. A coarse-grained mean-field approximation that resulted in a deterministic continuous-time flow model was introduced by [87], and its dynamics have been studied further [87,88].
Fig 11 illustrates the model. An mRNA consists of codons that are grouped into n sites, each site i has an associated occupancy level x i (t) 2 [0, 1] which can be interpreted as the probability that the site is occupied at time t. The ribosomes' inflow to the first site is λ 0 , which is known as the initiation rate, λ i is the elongation rate from site i to site i + 1, and λ n is the production rate. All rates are assumed to be positive. The ODE is written as follows:
The dynamics of the system above have been analyzed and shown to be monotone in [88]. In what follows, we provide an alternative approach that provides a Lyapunov function and establishes more powerful properties. Let y i ≔ 1 -x i , i = 1, .., n. Then, we can define a reaction network with species X i , Y i , i = 1, .., n as follows:
The network has 2n species, n + 1 reactions, and n conservation laws. It is depicted in Fig 11(b). The ODE system above describes the time-evolution of the reaction network with Mass-Action kinetics.
The graphical conditions of Theorem 4 are satisfied. Hence, ( 7) is an RLF for any n with R ¼ fR 1 ; R 2 ; :::; R nþ1 g. Since the network is conservative it follows that there exists a unique globally asymptotically stable steady state. Note that this results holds with general monotone kinetics.
In this subsection we show that our RLFs can provide valuable quantitative information regarding the behavior of the network beyond mere qualitative long-term behavior information.
Safety sets. Since our techniques are based on the construction of RLFs, we can compute safety sets which are the level sets of a Lyapunov function. If a system starts in a safety set it cannot leave it at any future time. Substituting Mass-Action kinetics, the safety set for a Lyapunov function Ṽ ðRðxÞÞ consists of piecewise polynomial surfaces and it is not necessarily convex. The safety set provided by an RLF surrounds all the steady states, i.e is not restricted to stoichiometric classes. In comparison, a concentration-dependent RLF provides a convex polyhedral safety set in a specific stoichiometric class. In order to illustrate this, consider the full PTM cycle with Mass-Action kinetics and let all the kinetic constants be 1. There are three conserved quantities, which we assume are set to Both safety sets corresponding to the two Lyapunov functions are chosen so that S = 2.5 lies on the boundary of the set. In other words, the substrate concentration is guaranteed not to exceed 2.5 if the system is initialized in the set. It can be clearly seen that the two sets are distinct, and they give different guarantees. Their intersection gives a tighter safety set.
Another example is a double processive PTM (Fig 8b) which has four dimensional stoichiometric classes. Hence, the 4D safety sets cannot be plotted, but their sublevel sets can still be visualized. Fig 12c and12d shows sublevel sets for different concentrations for the double phosphorylated species X 2 with total kinase, phosphatase and substrate concentrations fixed to 10AU each. Fig 12c ) shows the safety set with the concentration of the free kinase E not exceeding 2.5 and with [X 2 ] = 0. However, the sublevel set changes drastically if the concentration of X 2 is 0.5AU as shown in Fig 12d.
Flux analysis for the McKeithan network. Since the RLF are written in terms of rates (also called fluxes), our functions can be used in the context of flux analysis. Such techniques usually operate at steady state and do not take dynamics into consideration [89]. We provide an illustrative example to show how our RLF can be used. Let N = 2 for the network above. Usually, the network is initialized with zero concentration of the intermediate complexes.
The last inequality is included since the network is conservative and R 6 (m, ℓ) � R 6 ([M] T , [L] T ) holds due to the monotonicity of R.
The optimization problem above does not require knowledge of the kinetics as it is defined for fluxes. For the T-cell network, the solution of the problem is r � 6 ¼ 3r � 1 . This means that the flux r 6 is guaranteed to be less than 3r � 1 for all time. Converting these bounds to concentrations requires usage of the kinetics. Let R 1 (m, ℓ) = k 1 mℓ, and let R 6 ðc 2 Þ ¼ ac 2 1þbc 2 (Michaelis-Menten kinetics). Solving for c 2 , we can plot an upper bound on total amount of k 1 [M] T [L] T versus the maximum allowed concentration c 2 . If R 6 is Mass-Action then the relationship will be linear. Both curves are plotted in Fig 13.
We have presented a comprehensive theoretical framework and provided computational tools for the identification of a class of "structurally attractive" networks. It has been demonstrated that this class is ubiquitous in systems biology. Networks in this class have universal energylike functions called Robust Lyapunov Functions and, under additional mild conditions, can only admit unique globally stable steady states. Their Jacobians are well behaved and they cannot exhibit chaos, oscillations or multistability. The latter cannot be admitted even under inflow/outflow perturbations. Hence, LEARN can be used to rule out these networks as viable models for mechanisms that display such behaviors experimentally. Thus, our work supplements other mathematical methods used to invalidate models, as for example those in [90] and [91].
Our class of networks is distinct from the one identified by the HJF theory [14,17] and it has wider applications to biology as we have shown. Furthermore, our results include all networks that have been studied via compartmental system techniques [22,23] and via monotonicity techniques [8,24,33]. In fact, showing that the latter class of network always admits an RLF is a subject of a forthcoming paper. Refer to Table 1 for a comparison with techniques in the literature. In addition to wider applicability, our analysis has the advantage of showing persistence automatically, rather than needing to check it a priori as in [24]. Also, it has the advantage of having an explicit expression for the Lyapunov function which can be used for a deeper study of the dynamics such as the construction of safety sets and flux analysis as discussed before. In addition, Lyapunov functions have been extensively used to study the effect of interconnections, uncertainties, disturbances, and delays [9,10].
Our study of biochemical networks is not meant to be exhaustive, since we only focused on common motifs and cascades. We provide a computational package to help the wider community apply our techniques to study new networks.
We have presented the RLFs with two representations: rate-and concentration-dependent, and we have provided a toy example for dynamic flux analysis via a rate-dependent RLF. We look forward to these results being developed further to complement standard flux analysis techniques. For a given network, we have presented sufficient conditions for the existence of an RLF, and several necessary conditions. However, there are important networks that lie in the gap between the necessary and sufficient conditions. A relevant example is a ligand (L) binding a receptor (R), and initiating a PTM cycle for a substrate (S). The reaction network is:
It satisfies all necessary conditions but its global stability is still open. Future work includes the development of more general techniques to identify classes of networks that can be multi-stable but cannot admit oscillations or chaos. Furthermore, networks that admit RLFs have other strong properties in terms of contraction and stabilization [92], which will be studied in forthcoming papers.
We follow the standard notation and terminology on reaction networks [17,18,52,93].
A Chemical Reaction Network (CRN) consists of species and reactions. A species is what participates or is produced in a chemical interaction. In the context of biochemical networks a species can be a gene's promoter configuration, a substrate, an intermediate complex, an enzyme, etc. We denote the set of species by S ¼ fX 1 ; ::; X n g. A reaction is the transformation of reactants into products. Examples include binding/unbinding, decay, complex formation, etc. We denote the set of reactions by R ¼ fR 1 ; :::; R n g. Reactions have two distinct elements: the stoichiometry and the kinetics.
Stoichiometry. The relative number of molecules of reactants and products between the sides of each reaction is the stoichiometry. Hence, each reaction is customarily written as follows:
where α ij , β ij are nonnegative integers called stoichiometry coefficients. The expression on the left-hand side is called the reactant complex, while the one on the right-hand side is called the product complex. If a transformation is allowed to occur also in the opposite direction, the reaction is said to be reversible and its reverse is listed as a separate reaction. For convenience, the reverse reaction of R j is denoted as R -j . The reactant or the product complex can be empty, though not simultaneously. An empty complex is denoted by 0. This is used to model external inflows and outflows. An autocatalytic reaction is one which has a species appearing on both sides of the reaction simultaneously (e.g., D ! D + M). A network is called non-autocatalytic if it has no autocatalytic reactions.
The stoichiometry of a network can be summarized by arranging the coefficients in an augmented matrix n × 2ν as:
T called the stoichiometry matrix, which is defined as Γ = B -A, or element-wise as:
Kinetics. The relations that determine the velocity of transformation of reactants into products are known as kinetics. We assume an isothermal well-stirred reaction medium. In order to study kinetics, a nonnegative number x i is associated to each species X i to denote its concentration. Assume that the chemical reaction R j takes place continuously in time. A reaction rate or velocity function We do not assume particular kinetics. We only assume that the reaction rate functions R j (x), j = 1, ..ν satisfy the following minimal assumptions: AK1. each reaction varies smoothly with respects to its reactants, i.e R j (x) is continuously differentiable;
AK2. each reaction needs all its reactants to occur, i.e., if α ij > 0, then x i = 0 implies R j (x) = 0;
AK3. each reaction rate is monotone with respect to its reactants, i.e @R j /@x i (x) � 0 if α ij > 0 and @R j /@x i (x) � 0 if α ij = 0;
AK4. The inequality in AK3 holds strictly for all positive concentrations, i.e when x 2 R n þ . Reaction rate functions satisfying AK1-AK4 are called admissible. For given stoichiometric matrices A, B, the set of admissible reactions is denoted by K A .
Dynamics. The dynamics have been already given in (1). The set
þ is forward invariant for any initial condition x � , and it is called the stoichiometric compatibility class associated with x � . For a conservative network all stoichiometric classes are compact convex polyhedral sets.
We sometimes will use the following assumption which is necessary for the existence of positive steady states.
AS1. There exists v 2 ker Γ such that v � 0.
We have provided an informal definition of the notion of RLF in the introduction. The inequality in Eq (4) must hold for all R 2 K A . As observed before, AK1-AK4 imply a zero-sign pattern on @R/@x (see Fig 2d for an illustration). This motivates defining the class of matrices with the specific sign pattern as follows:
where P is the set of reaction-reactant pairs defined before. Definition 1. Given a network ðS ; RÞ. A locally Lipschitz function Ṽ : R n ! R �0 is said to be an RLF if it satisfies the following:
2. D Ṽ :¼ ð@ Ṽ =@rÞrGr � 0 for all r 2 K A and all r for which @ Ṽ =@rðrÞ exists.
At points of non-differentiability, the time-derivative of VðxÞ ¼ Ṽ ðRðxÞÞ is defined in the sense of Dini (see S1 Text §1.1 for a review of Lyapunov theory and generalized derivatives).
We will show how the rank-one matrices Q 1 , .., Q S (defined in the Results section) can be used to embed the dynamics of the nonlinear network in a cone of linear systems. Although the Lyapunov function Ṽ ðRðxÞÞ is a function in the concentration x, it is defined as a composition V ¼ Ṽ � R. Therefore, we study the ODE in reaction coordinates. Let x(t) be a trajectory that satisfies (1) and let r(t)≔ R(x(t)). Hence,
where rðtÞ ≔ @R @x ðxðtÞÞ 2 K A . The basic idea is to consider ρ(t) as an unknown time-varying matrix. Since its zero-sign pattern is known, we can decompose ρ(t) in the following way:
where [ρ(t)] ji = ρ ji (t) > 0, and [E ji ] j 0 i 0 = 1 if (j 0 , i 0 ) = (j, i) and zero otherwise. The matrices {E ji |(i, j) such that α ij = 0} form the canonical basis of the matrix space K A . Substituting (18) in (17) we can embed the dynamics of the network (17) in the conic combinations of a finite set of extremal linear systems as follows:
where Q ℓ , ℓ = 1, .., s have been defined before Theorem 1. This also implies that the Jacobian of ( 17) can be written at any interior point as: ð@R=@xÞG ¼ P s '¼1 r ' Q ' . Hence, the Jacobian belongs to cone generated by the extremals Q 1 , .., Q s . Note that (19) can be interpreted as representing a linear parameter-varying system which has s nonnegative time-varying parameters {ρ 1 (t), .., ρ s (t)}. The linear systems are given by rank-one extremals Q 1 , .., Q s . The proof of Theorem 1 is completed in S1 Text §1.2.
The results presented in [26,27] have been derived via a direct analysis of the associated reaction networks. The framework introduced above enables interpreting these results in a more general framework and allows generalizing them. Hence we revisit the algorithms introduced for the existence and construction of PWL RLFs, and implement them in the LEARN MATLAB package. Furthermore, we also introduce piecewise quadratic RLFs based on the new framework introduced in this paper.
Piecewise linear RLFs. Consider a CRN (1) with a G 2 R n�r and a given partitioning matrix H 2 R p�r such that ker H = ker Γ. A PWL RLF is piecewise linear-in-rates, i.e., it has the form: VðxÞ ¼ Ṽ ðRðxÞÞ, where Ṽ : R n ! R is a continuous PWL function. Assuming AS1, the piecewise linear function is given as Construction via linear programming. Based on Theorem 1, we present a simpler linear program than the one presented in [27]. The proof is presented in S1 Text §2.2.
Theorem 2. Given a network ðS ; RÞ that satisfies AS1 and a partitioning matrix H 2 R p�r . Let {v i } be a basis for ker Γ. Consider the linear program:
Then there exists a PWL RLF with partitioning matrix H if and only if there exists a feasible solution to the above linear program that satisfies ker C = ker Γ.
Remark 1. The linear program above does not enforce convexity on Ṽ . Nevertheless, LEARN allows the user to search amongst convex Ṽ 's only. See S1 Text §2. 3.
In LEARN there is a default choice for the matrix H, and it also allows for a manual input by the user. The default choice is H = Γ which gives the following Lyapunov function (where the SoC RLF introduced in ( 6) is a special case):
The user can add rows to H. Usually rows of the form {γ i ± γ j |i, j = 1, .., n}, where γ 1 , .., γ n are the rows of Γ, are good candidates.
Networks without positive steady states. If AS1 is not satisfied, then a linear program can be designed for constructing RLFs over a given partition. This is discussed in S1 Text §2.4.
An iteration for the construction of convex PWL RLFs. Assuming both AS1 and allowing non-autocatalytic networks only, a computationally-light iterative algorithm for constructing a convex Lyapunov function was presented in [26,27]. Here we generalize the algorithm by dropping these two assumptions. The objective is to find a matrix C ¼ ½c T 1 ; ::::; c T m � T such that Ṽ ðrÞ ¼ max k¼0;::;m c T k r is a Lyapunov function, where c 0 ≔ 0. We state the algorithm below. We use the notation supp(c k ) = {R j |c kj 6 ¼ 0}, which is the set of all those reactions that appear in c T k r, and let IðR j Þ ¼ fX i ja ij > 0g which is the set of reactants for reaction R j . We have the following result, which is proved in S1 Text §2.5.
Theorem 3. Given a network ðS ; RÞ. Let G ¼ ½g T 1 ; :::; g T n � T 2 R n�n be its stoichiometry matrix. If the following algorithm terminates successfully, then Ṽ is an RLF.
Parameters: N as the upper maximum number of iterations.
Success. Ṽ ðrÞ ¼ max k¼0;::;m c T k r is the desired function else
The algorithm did not converge within the prescribed upper maximum number of iterations. end
The algorithm above is computationally very light compared to the linear program with a large H. Furthermore, if the network satisfies AS1 then the RLF can be written as Ṽ ðrÞ ¼ k Cr k 1 .
Graphical criteria for the construction of Max-Min RLFs. Compared to computational conditions, it is highly desirable to have graphical conditions and some have been provided in [26,27]. We reformulate those conditions to be more friendly for computational implementation in LEARN. Those conditions enable the identification of attractive networks by mere inspection of the reaction graph for a particular class of networks.
We introduce some notations. Let ðS ; RÞ be a given non-autocatalytic network that satisfies AS1. Consider the decomposition R ¼ R r [ R i into the subsets of reactions that are reversible and irreversible, respectively. Furthermore, we can decompose R
r Þ be the corresponding irreversible subnetwork and let~be its stoichiometry matrix. Since the designation of a forward and reverse reaction is arbitrary, we need a decomposition such that G has a one-dimensional nullspace. If such a decomposition exists, then we call the original network ðS ; RÞ an M-network. Our graphical condition applies to this class of networks, and it can be stated as follows.
Theorem 4. Let ðS ; RÞ be an M-network, and let ðS ; R þ r [ R i Þ be the subnetwork defined above, where the reactions are enumerated as R þ r ¼ fR 1 ; :::; R n 1 g, R i ¼ fR n 1 þ1 ; :::; R ñ g. If the irreversible subnetwork satisfies the following properties:
1. each species participates in exactly one reaction, and
r satisfies the following statement: If a species X i is a product of R j , then X i is not a product of another reaction, then Ṽ ðRðxÞÞ ¼ max RðxÞ min RðxÞ;
, is a convex PWL RLF, where w = [w 1 , . . ., w |ν 1 |] T belongs to the null space of G. Piecewise quadratic-in-rates RLFs. The framework developed in this paper allows us to go beyond PWL RLFs, and consider other classes of functions such piecewise quadratic-inrate functions of the form:
for some matrices
Instead of linear programming, construction of PWQR RLFs is a copositive programming problem. Although copositive programs are convex, solving them generally is shown to be NPhard [94]. Therefore, we use a common relaxation scheme based on the observation that the class of copositive matrices encompasses the classes of positive semi-definite matrices, and nonnegative matrices. The following theorem states the result and it is proven in S1 Text §3.1.
Theorem 5. Given a network ðS ; RÞ that satisfies AS1 and a partitioning matrix H 2 R p�r . Let {v i } be the basis for the kernel of Γ Consider the following semi-definite program: Find
Robust non-degeneracy. It has been shown in [27] that the negative Jacobian of any network admitting a PWL RLF is P 0 , which means that all principal minors are nonnegative. We show that the reduced Jacobian (i.e., Jacobian with respect to a stoichiometric class) is nondegenerate for all admissible kinetics if it is so at one interior point only. The proof is stated in S1 Text §4.1.
Theorem 7. Assume that there exists a PWL RLF. If for some kinetics R 2 K A there exists a point in the interior of a proper stoichiometric class such that the reduced Jacobian is non-singular at it, then the reduced Jacobian is non-singular in the interior of R n þ for all admissible kinetics.
In LEARN, robust non-degeneracy is checked with ρ ℓ = 1, ℓ = 1, . . ., s. It amounts to checking the non-singularity of one matrix. [27], which is automatically satisfied for conservative M-networks. Alternatively, global stability follows automatically for any positive steady state if the network is robustly nondegenerate [95]. Hence, Theorem 7 can be used to verify global stability when a PWL RLF exists. Note, however, that the test above is with respect to the stoichiometric class only. In the case of degenerate reduced Jacobians, a stoichiometric class can be partitioned further into kinetic compatibility classes [16]. The graphical LaSalle's algorithm applies to such cases also.
Absence of critical siphons. A siphon is any (minimal) set of species which has the following property: if those species start at zero concentration, then they stay so during the course of the reaction [41]. Siphons are of two types: trivial and critical. A trivial siphon is a siphon that contains the support of a conservation law. A critical siphon is a siphon which is not trivial. Critical siphons can be found easily from the network graph. The absence of critical siphons in a network has been shown to imply that it is structurally persistent (for conservative networks or systems with bounded flows) [41]. Informally, a system is persistent if the following holds: if all species are initialized at nonzero concentrations, none of them will become asymptotically extinct. We show that the existence of critical siphons precludes the existence of RLF under mild conditions which serves as an easy-to-check condition to preclude the existence of an RLF. Review of the concept of siphons and the proof the result is included in S1 Text §4.2.
Theorem 8. Given a network ðS ; RÞ that satisfies AS1. Assume it has a critical siphon P � S . Let LðPÞ � R be the set of reactions for which the species in P are reactants. Then there cannot exist a PWL RLF if any of the following holds:
1. lðPÞ ¼ R, i.e P is a critical deadlock.
In this subsection we study an alternative RLF and we link the results with the ones proposed in [34,36]. We will show that any RLF has an alternative form if it satisfies a mild condition.
In particular, all PWL RLFs have alternative forms. Assume that (1) has a steady state x e . Then, we ask whether there exists a Lyapunov function of the form VðxÞ ¼ V ðx x e Þ. However, note that this Lyapunov function decreases only in the stoichiometric class containing x e and that computing its level sets requires knowing x e . We call V a concentration-dependent RLF.
Similar to before, we will characterize the existence of an RLF of the form V ðx x e Þ for a network ðS ; RÞ by the existence of a common Lyapunov function for a set of extremals of an appropriate cone. In this subsection, we assume that there exists a positive r 2 R n þ such that Γr = 0.
We will adopt an alternative representation of the system dynamics. Consider a CRN as in (1), and let x e be a steady state. Then, there exists x@ðxÞ 2 R n þ such that (1) can written equivalently as:
The existence of x@ ≔ x e + ε x (x -x e ) for some ε x [0, 1] follows by applying the Mean-Value Theorem to R(x) along the segment joining x e and x.
Similar to the analysis for a rate-dependent RLFs, the Jacobian of (1) can be shown to belong to the conic span of a set of rank-one matrices fG i 1 e T j 1 ; :::; G i s e T j s g where {Γ 1 , .., Γ n } are the columns of Γ. The pairs (i ℓ , j ℓ ), ℓ = 1, .., s are the same pairs used before.
Let D T be a matrix with columns that are the basis vectors of ker Γ T . The following theorem is proven in S1 Text §5.1.
Theorem 9. Given a network ðS ; RÞ. There exists a common Lyapunov function V : R n ! � R þ for the set of linear systems f _ z ¼ ðG i 1 e T j 1 Þz; :::; _ z ¼ ðG i s e T j s Þzg, on the invariant subspace {z: D T z = 0} if and only if V ðx x e Þ is a concentration-dependent RLF for any x e .
Relationship between the RLFs in concentration and rates. We show next that if Ṽ is a rate-dependent RLF that satisfies a relatively mild additional assumption, then then V ðx x e Þ is a concentration-dependent RLF, where x e is a steady state point for (1). The following theorem can be stated and is proved in S1 Text §5.2.
Theorem 10. Let Ṽ be an RLF for the network ðS ; RÞ. If there exists V : R n ! � R þ such that for all r 2 R n :
then V is a concentration-dependent RLF for the same network. PWL functions in concentrations. All PWL RLFs constructed before have the property that there exists V such that Ṽ ðrÞ ¼ V ðGrÞ. Hence, there exists a concentration-dependent PWL RLF for the same network. In particular, consider a PWL RLF defined with a partitioning matrix H as in (20). By AS1 and the assumption that ker H = ker Γ, there exists G 2 R p�n and B 2 R where it can be seen that V k has nonempty interior iff W k has nonempty interior.
Therefore, as the pair (C, H) specifies a PWL RLF, the pair (B, G) also specifies the function:
V ðzÞ ¼ b T k z; when S k Gz � 0;
where B ¼ b 1 ; :::; bm 2 h i T . If Ṽ is convex, then it can be written in the form: V 1 (x) = kCR(x)k 1 .
Similarly, the convexity of V implies that V 2 (x) = kB(x -x e )k 1 , where the latter is the Lyapunov function used in [36].
Theorem 10 shows how to go from a rate-dependent to a concentration-dependent RLF.
PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1007681February 24, 2020