A computational framework for a Lyapunov-enabled analysis of biochemical reaction networks

Complex molecular biological processes such as transcription and translation, signal transduction, post-translational modification cascades, and metabolic pathways can be described in principle by biochemical reactions that explicitly take into account the sophisticated network of chemical interactions regulating cell life. The ability to deduce the possible qualitative behaviors of such networks from a set of reactions is a central objective and an ongoing challenge in the field of systems biology. Unfortunately, the construction of complete mathematical models is often hindered by a pervasive problem: despite the wealth of qualitative graphical knowledge about network interactions, the form of the governing nonlinearities and/or the values of kinetic constants are hard to uncover experimentally. The kinetics can also change with environmental variations. This work addresses the following question: given a set of reactions and without assuming a particular form for the kinetics, what can we say about the asymptotic behavior of the network? Specifically, it introduces a class of networks that are “structurally (mono) attractive” meaning that they are incapable of exhibiting multiple steady states, oscillation, or chaos by virtue of their reaction graphs. These networks are characterized by the existence of a universal energy-like function called a Robust Lyapunov function (RLF). To find such functions, a finite set of rank-one linear systems is introduced, which form the extremals of a linear convex cone. The problem is then reduced to that of finding a common Lyapunov function for this set of extremals. Based on this characterization, a computational package, Lyapunov-Enabled Analysis of Reaction Networks (LEARN), is provided that constructs such functions or rules out their existence. An extensive study of biochemical networks demonstrates that LEARN offers a new unified framework. Basic motifs, three-body binding, and genetic networks are studied first. The work then focuses on cellular signalling networks including various post-translational modification cascades, phosphotransfer and phosphorelay networks, T-cell kinetic proofreading, and ERK signalling. The Ribosome Flow Model is also studied.


Introduction
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].

Overview and comparison
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.

Terminology and motivational example
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 1 ; :::; R n g. (see Methods for an elaborate discussion) Fig 2a) gives an example of a reaction network for a core signaling motif which is the standard post-translational modification (PTM) cycle [48,49]. The relative gain or loss of molecules of each species in a reaction is encoded in a matrix G 2 R n�n called the stoichiometry matrix. It is given in Fig 2b for the PTM cycle. CRNs admit graphical representations naturally. A CRN can be modeled as a graph with two sets of nodes: reactions and species. Mathematically, it is a bipartite weighted directed graph, called the species-reaction graph (or a Petri-net [50]). The graph corresponding to the Table 1. 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]  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.

Robust Lyapunov functions
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: for all x and for all R admissible: 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).

Characterization of RLFs
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, ., γ 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

The search for RLFs
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: Two special cases are of interest to us: Sum-of-currents (SoC) RLFs. These are functions of the form: where x ¼ ½x 1 ; ::; x n � 2 R n �0 is a positive vector and k[z 1 , .., z n ] T k 1 ≔ ∑ i |z i | is the 1-norm. The function considered in [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 functionV such thatV ðx À x e Þ (where x e is a steady state) is a Lyapunov function for any admissible R. We callV a concentration-dependent RLF. We show thatṼ ðrÞ ¼ k BGr k 1 is an RLF iffV ðzÞ ¼ k Bz k 1 is a concentration-dependent RLF (see Methods, Theorem 11). These PWL functions relate to the ones proposed in [34,36]. Note, however, thatV ð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 class of structurally attractive biochemical networks
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.

Binding/Unbinding reactions
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.
Simple binding reaction. Fig 3a represents a simple reversible binding reaction: 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. Simple binding with enzyme inflow-outflow. Fig 3b represents the following binding reaction with enzyme inflow-outflow: 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:

XE n
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 ! E, 0 ! X, nE ! E n , E n + X ! XE n , XE n ! 0. Competitive binding reaction. The following reactions (depicted in Fig 3d) describe the situation when two molecules E 1 , E 2 compete to bind with X: 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

Three-body binding
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 networks
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 Rib þ mRNA Ð mRNA : Ribo ! mRNA þ Ribo; mRNA : Ribo ! Rib and 0 ! P ! 0. Applying Theorem 3 to the first network we get the following Lyapunov function:Ṽ ðRðxÞÞ ¼ max fR 4 ðxÞ; R 1 ðxÞ À R 2 ðxÞ À R 3 ðxÞ À R 4 ðxÞ; À R 1 ðxÞ þ R 2 ðxÞ þ R 3 ðxÞg: 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 enzymatic networks
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: VðxÞ ¼ jR 1 À R À 1 ðxÞ À R 2 ðxÞj þ jR 2 ðxÞ À R 3 ðxÞj þ jR 1 ðxÞ À R À 1 ðxÞ À R 3 ðxÞj; 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: A À À ! R 6 A; P þ À ! R 7 P, 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.

Post-translational modification cycle cascades
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]. Fig 8b) depicts a multiple-site futile cycle with a processive mechanism. 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 R ¼ fR k À R À k ; k ¼ 1; ::; ng, and R −k (x) :� 0 if R k is irreversible.
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. (Similar observations have been made in [72][73][74].) Fig 1 depicts sample trajectories for the processive and distributive cycle with Mass-Action kinetics.

Phosphotransfer and phosphorelay networks
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 take 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 R ¼ fR 1 À R À 1 ; R 2 À R À 2 ; R 3 À R À 3 ; R 4 ; R 5 À R À 5 ; R 6 g. A phosphorelay. 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.

T-cell kinetic proofreading network
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 ,  [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.  [85]. It can be described using the network:

ERK signaling pathway with RKIP regulation
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.

The ribosome flow model
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.

Quantitative analysis via RLFs
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   (a),(b), Safety sets for the PTM cycle (Fig 6c). (a) The safety set corresponding to the rate-dependent RLF for the PTM cycle. It is the α-level set of V where α has been chosen such that the concentration of S does not exceed 2.5. (b) The safety set corresponding to the concentration-dependent RLF. The safety set has been chosen similarly to satisfy the same condition. (c),(d), Sub-levels sets for the safety sets corresponding to the rate-dependent RLF (7) for the double processive PTM cycle (Fig 8b).  Fig 12b. 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 and 12d 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.

Discussion
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.
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: The two submatrices A, B can be subtracted to yield an n × ν matrix G ¼ ½g T 1 ::g T n � 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 R j : � R n þ ! � R þ is assigned to each reaction. The widely-used Mass-Action kinetics have the following expression: R j ðxÞ ¼ k j Q n i¼1 x a ij i , where k j , j = 1, .., ν are positive numbers known as the kinetic constants. Many other kinetic forms are used in biology such as Michaelis-Menten, Hill kinetics, etc.
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 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 C x� ≔ ðfx� g þ ImðGÞÞ \ � R n þ 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.

RLFs and the decomposition of the dynamics
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.

Computational construction of RLFs
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 where the regions W k ¼ fr 2 R n : S k Hr � 0g; k ¼ 1; ::; m form a proper conic partition of R n , while fS k g Verifying a candidate RLF. Checking if a given PWL function is an RLF can be posed as a linear program. It is discussed in S1 Text §2.1 and is coded into LEARN.
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: 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. 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 ¼ R þ r [ R À r into the forward and backward reactions, respectively. Let ðS ; R i [ 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 thatG 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 2. each reaction R j 2 R þ 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Þ; where R ¼ 1 w 1 ðR 1 À R À 1 Þ; :::; 1 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 P k 2 R n�n ; c k 2 R n , k = 1, .., m.
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 where d = dim(ker Γ) and N k is the set of neighbor of region W k (see S1 Text §3.2). If the SDP is feasible, thenṼ as defined in (22) is an RLF for ðS ; RÞ if kerṼ ¼ ker G. This class of networks for which PWQ RLFs exist is potentially larger than that of PWL RLFs even when we set c k = 0, k = 1, .., m in (22) as the following proposition establishes. The proof is given in S1 Text §X.

Properties of attractive networks
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.

Remark 2.
Robust non-degeneracy, coupled with the existence of a PWL RLF, automatically guarantees the uniqueness of positive steady states and their exponential stability (see S1 Text §4. 1.2, §4.1.3). Globally stability has been checked via a LaSalle algorithm in [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.

ðS ; RÞ is a conservative M network.
3. ðS ; RÞ is conservative and has a positive non-degenerate steady state for some admissible kinetics.

RLFs in other coordinates
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 callV a concentration-dependent RLF. Similar to before, we will characterize the existence of an RLF of the formV ð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 2 [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 ifV ð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 thenV ð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 existsV : R n ! � R þ such that for all r 2 R n :Ṽ

ðrÞ ¼V ðGrÞ; ð28Þ
thenV is a concentration-dependent RLF for the same network. PWL functions in concentrations. All PWL RLFs constructed before have the property that there existsV 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 ofV 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. The following theorem shows that one can start with either PWL RLF to get the other. It is proved in S1 Text §5.3.

Remark 4.
Since D T (x − x e ) = 0 for x 2 C x e , then if kB(x − x e )k 1 is an RLF, then k(B + YD T ) (x − x e )k 1 is also an RLF for an arbitrary matrix Y. Furthermore, since Theorem 11 has shown that the concentration-based and the rate-based representations are equivalent, it is easier to check and construct RLFs in the rate-based formulation and they hold the advantage of being decreasing for all trajectories over all stoichiometry classes.
Computational package. Calculations were performed using MATLAB 10 via our software package LEARN available at https://github.com/malirdwi/LEARN. Available subroutines and example runs are included in S1 Text §7. The package cvx [96] has been used for solving linear and semi-definite programs, and the package PetriBaR for enumerating siphons [97].
Supporting information S1 Text. Supporting information file with mathematical proofs, generalization of the results and additional information. (PDF)