Robust Lyapunov Functions for Reaction Networks: An Uncertain System Framework
aa r X i v : . [ m a t h . O C ] O c t Robust Lyapunov Functions for Reaction Networks:An Uncertain System Framework
M. Ali Al-Radhawi a , David Angeli b,c a Department of Mechanical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge,MA 02139, United States. b Dept. of Electrical and Electronic Engineering, Imperial College London, London SW7 2AZ, United Kingdom. c Dip. di Ingegneria dell’Informazione, University of Florence, 50139 Florence, Italy.
Abstract
We present a framework to transform the problem of finding a Lyapunov function for a ChemicalReaction Network (CRN) with arbitrary monotone kinetics expressed in concentration coordi-nates into finding a common Lyapunov function for a linear differential inclusion in reactioncoordinates. Alternative formulations in different coordinates are also provided. This is appliedto reinterpret previous results by the authors on Piecewise Linear in Rates Lyapunov functionsand to establish a link with contraction analysis. The framework is then applied to derive power-ful results on the network’s persistence and the uniqueness of equilibria.
Keywords:
Lyapunov stability, Persistence, Contraction Analysis, Linear DifferentialInclusions, Biochemical Networks.
1. Introduction
Chemical Reaction Networks (CRNs) are used as models in several areas of science and engi-neering; including chemical engineering, population dynamics, and molecular systems biology.The last application has drawn the recent interest of the systems and control community to thisfield [1, 2].One of the central dilemmas in systems biology is that the kinetic information required toconstruct a detailed mathematical model is scarce and subject to uncertainty. This is unlike thewealth of information available on the graphical description of the networks involved. Hence, themain focus of Chemical Reaction Networks Theory and of our research has been to draw con-clusions about the asymptotic behaviour of networks classes regardless of detailed knowledge ofits kinetics. This has been advocated towards the aim of constructing “complex biology without
Email addresses: [email protected] (M. Ali Al-Radhawi), [email protected] (DavidAngeli)
Preprint submitted to Elsevier February 27, 2018 arameters” [3]. The presumed feasibility of this goal has been motivated by the observationthat large classes of CRNs converge asymptotically to steady states regardless of the model ofkinetics involved. Indeed, partial success has been achieved in this pursuit. For example, theclass of weakly reversible zero-deficiency networks with Mass-Action kinetics has been shownto have unique asymptotically stable equilibria in the interior of the orthant [4, 5]. Monomolecu-lar networks, on the other hand, have also been handled within the framework of compartmentalsystems [6]. The theory of monotone systems has been also applied to provide graphical condi-tions for global convergence of some of these networks [7]. Recently, piecewise linear Lyapunovfunctions based methods has been proposed recently [8, 9, 10, 11, 12].In a previous work [10, 11], the authors proposed a direct approach to the problem, wherePiecewise-Linear in Rates (PWLR) Lyapunov functions have been introduced. In addition totheir simple structure, these functions are robust with respect to arbitrary variations of kineticconstants, and only require mild assumptions on the reaction rates, with mass-action kineticsbeing a special case.In this paper, we generalize this approach by transforming the problem of finding a Lya-punov function expressed in concentration coordinates for networks with arbitrary monotonekinetics into finding a common Lyapunov function for a linear parameter varying system writtenin reaction coordinates. Furthermore, several related Lyapunov functions are presented. As aresult, we link the PWLR Lyapunov functions introduced in [10] with results known in literaturefor piecewise linear Lyapunov functions [13, 14, 15]. Furthermore, we interpret our results interms of contraction analysis and variational dynamics. The existence of Lyapunov function hasalso strong algebraic and dynamical implications. Particularly, results on the persistence and theuniqueness of equilibria are also established.This paper is organized as follows. In Section 2, we present the background and assump-tions. Section 3 defines what we mean by a Robust Lyapunov function, and develops the uncer-tain systems framework. PWLR Lyapunov functions are introduced in Section 4, uniqueness ofequilibria and persistence are presented in Section 5. Section 6 discusses the relationship withcontraction analysis. Proofs are collected in the Appendix.Note that some of the results in Section 3 and 4 in this paper have been published preliminarilywithout proofs in [16].
Notation.
Let A ⊂ R n be a set, then A ◦ , ¯ A, ∂A, co A denote its interior, closure, boundary,and convex hull, respectively. Given x ∈ R n , a column vector, its ℓ ∞ -norm is k x k ∞ =max ≤ i ≤ n | x i | . The inequalities x ≥ , x > , x ≫ denote elementwise nonnegativity,elementwise nonnegativity with at least one positive element, and elementwise positivity, re-spectively. A is a signature matrix if it is diagonal and the diagonal entries belong to {± } . For A ∈ R n × ν , ker( A ) denotes the kernel or null-space of A , while Im ( A ) denotes the image spaceof A . A ∈ R n × n is Metzler if all off-diagonal elements are nonnegative. The set of n × n S n . Let A ∈ S n , then A ≥ ( > )0 denotes A being posi-tive semi-definite (definite), respectively. A (cid:23) denotes entrywise nonnegativity. The all-onesvector is denoted by , where its dimension can normally be inferred from the context. Let { A i } ki =1 ⊂ R n × m , its conic hull denotes the set { P ki =1 λ i A i : λ i ∈ ¯ R + } . Let V : D → R , thenthe kernel of V is ker( V ) = V − (0) . T M is the tangent bundle of the differentiable manifold M .
2. Background on Chemical Reaction Networks
The field of CRN dynamics has an established literature [5, 2]. We review here the relevantbackground. A reaction network has two mathematical defining features: the stoichiometry andthe kinetics . Informally speaking, stoichiometry describes the relative number of molecules ofreactants and products involved whenever each reaction occurs, while kinetics is concerned withthe relations that govern the velocity of transformation of reactants into products. We explainboth below.
A Chemical Reaction Network (CRN) is defined by a set of species S = { X , .., X n } , anda set of reactions R = { R , ..., R ν } . Each reaction is denoted as: R j : n X i =1 α ij X i −→ n X i =1 β ij X i , j = 1 , .., ν, (1)where α ij , β ij are nonnegative integers called stoichiometry coefficients . The expression on theleft-hand side is called the reactant complex , while the one on the right-hand side is called the product complex. The forward arrow refers to the idea that the transformation of reactants intoproducts is only occurring in the direction of the arrow. In order to allow external inflows oroutflow the reactant or product complex can both be empty, though not simultaneously.The stoichiometry of a network can be summarized by arranging the coefficients in an aug-mented matrix n × ν as: ˜Γ = [ A | B ] , where [ A ] ij = α ij , [ B ] ij = β ij . (2)The two matrices can be subtracted to yield an n × ν matrix Γ = [ γ T ..γ Tn ] T called the stoichiometry matrix , which is defined as Γ = B − A , or element-wise as: [Γ] ij = β ij − α ij . A left null vector d ∈ R n , d T Γ = 0 with d > is said to be a conservation law . If thereexists a conservation law d ≫ , the network is said to be conservative .3 reaction R j is said to be an input reaction to X i if β ij > , and is said to be an output reaction of X i if α ij > . Let P ∈ S be a non-empty set of species. Denote the set of output reactionsof species in P by Λ( P ) . Then, a nonempty set P ⊂ S is called a siphon [17] if each inputreaction to a species in P is also an output reaction of species in P . A siphon is a deadlock if Λ( P ) = R . A siphon or a deadlock is said to be critical if it does not contain a set of speciescorresponding to the support of a conservation law. Assume we have an isothermal well-stirred chemical reactor; this implies that the speciesare distributed uniformly in the reactor. In order to study kinetics, a nonnegative number x i isassociated 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 assignedto each reaction.A widely-used expression which originates from statistical thermodynamics is given as: R j ( x ) = k j n Y i =1 x α ij i , (3)(the so called Mass-Action kinetics), with the convention = 1 , where k j , j = 1 , .., m arepositive numbers known as the kinetic constants, and are usually highly uncertain.The Mass-Action model is not a universal model since there are other models which are popularin systems biology like the Michaelis-Menten model, and the Hill model.In this work we will pursue ”kinetics-independent” approach. More precisely, we assume thatthe reaction rate function of a CRN is unknown except for satisfying the following assumptions: AK1. it is a C function, i.e. continuously differentiable; AK2. x i = 0 ⇒ R j ( x ) = 0 , for all i and j such that α ij > ; AK3. it is nondecreasing with respect to its reactants, i.e ∂R j ∂x i ( x ) ( ≥ α ij >
0= 0 : α ij = 0 . (4) AK4.
The inequality in (4) holds strictly for all x ∈ R n + .Reaction rate functions satisfying AK1-AK4 are called admissible . For a given stoichiomet-ric matrices A, B , the set of admissible reactions is denoted by K A . A network family is thetriple N A,B = ( S , R , K A ) . Remark 1.
The assumptions above imply the monotonic dependence of the reaction rate on theconcentration of its reactants. This captures the basic intuition about the nature of a reaction ince, as the concentration of reactants increases, the likelihood of collision between moleculesincreases, and hence the rate of the reaction. Note that, in principle, a monotonically decreasingdependence on the reactants can also be accommodated to model inhibition . Remark 2.
It can be noted that if a reaction rate function is admissible, then the correspondingJacobian ∂R/∂x exhibits a certain zero sign-pattern that can be read from the graph, and anymatrix satisfying that pattern correspond to an admissible R .2.3. Dynamics The dynamics of a CRN with n species and ν reactions are described by a system of ordinarydifferential equations (ODEs) as: ˙ x ( t ) = Γ R ( x ( t )) , x (0) ∈ ¯ R n + (5)where x ( t ) is the concentration vector evolving in the nonnegative orthant ¯ R n + , Γ ∈ R n × ν is thestoichiometry matrix, R ( x ( t )) = [ R ( x ( t )) , R ( x ( t )) , ..., R ν ( x ( t ))] T ∈ ¯ R ν + is the reaction ratesvector.Note that (5) belongs to the class of positive systems , i.e, ¯ R n + is forward invariant. In addition,the manifold C x ◦ := ( { x ◦ } + Im (Γ)) ∩ ¯ R n + is forward invariant, and it is called the stoichiometriccompatibility class associated with x ◦ . Therefore, all stability results in this paper are relativeto the stoichiometry compatibility class. Note that for a conservative network all stoichiometricclasses are compact convex polyhedral sets.An equilibrium x e of (5) is non-degenerate if the Jacobian evaluated at x e relative to C x e is nonsingular. More precisely, considering changing coordinates using a transformation matrix T = [ T T D ] T , where D T has full row rank and D T Γ = 0 , and T is any matrix such that T isnonsingular. Then, the Jacobian in the new coordinates can be written as: T Γ ∂R∂x T − = " J J . (6)Therefore, x e is nondegenerate iff J evaluated at x e is nonsingular. J is called a reducedJacobian .The stoichiometry of the network will be assumed to satisfy the following assumption: AS There exists v ∈ ker Γ such that v ≫ . This condition is necessary for the existence of asteady state in which all concentrations are positive.5 . Robust Lyapunov Functions and Linear Differential Inclusions In order for the stability analysis of CRNs to be independent of the specific kinetics, weaim at constructing Lyapunov functions which are dependent only on the graphical structure,and hence are valid for all reaction rate functions that belongs to K A . Therefore, we state thefollowing definition. Definition 1 (Robust Lyapunov Function) . Consider (5) and let x e be an equilibrium. Let ˜ V :¯ R q → ¯ R + be locally Lipschitz, and let W R,x e : R n → R q be a C function. Then, ( ˜ V , W
R,x e ) is said to induce a Robust Lyapunov Function (RLF) with respect to the network family N A,B iffor any choice of R ∈ K A , x e ∈ ¯ R n + , the function V R,x e = ˜ V ◦ W R,x e is
1. Positive-Definite : V R,x e ( x ) ≥ , and V R,x e ( x ) = 0 if and only if R ( x ) ∈ ker Γ .
2. Nonincreasing : ˙ V R,x e ( x ) ≤ for all x ∈ C x e .A network for which an RLF exists is termed a Graphically Stable Network (GSN).
Remark 3.
As will be seen later, the function ˜ V used in the definition of the Lyapunov function(through composition) is invariant with respect to the specific network realization in K A , whilethe function W R,x e is allowed to depend on the kinetics of the network. Two main examples ofthe function W R,x e are W R,x e ( x ) = R ( x ) , and W R,x e ( x ) = x − x e . With a mild abuse ofterminology, we call the parameterized Lyapunov function V R,x e an RLF. Remark 4.
The time-derivative in the definition above is the upper right Dini’s derivative [18]: ˙ V ( x ) := lim sup h → + V ( x + h Γ R ( x )) − V ( x ) h , (7) which is finite for all x since V is locally Lipschitz. Since the RLF defined above is not strict, we need the following definition.
Definition 2 (The LaSalle’s Condition) . An RLF V R,x e for N A,B is said to satisfy the LaSalle’sCondition if for any choice R ∈ K A the following statement holds.If a solution ϕ ( t ; x ◦ ) of (5) satisfies ϕ ( t ; x ◦ ) ∈ ker ˙ V ∩ C x e , t ≥ , then this implies that ϕ ( t ; x ◦ ) ∈ E x ◦ for all t ≥ , where E x ◦ ⊂ C x ◦ is the set of equilibria for (5) contained in C x ◦ . The following theorem adapts Lyapunov’s second method [18, 11] to our context.
Theorem 1 (Lyapunov’s Second Method) . Given (5) with initial condition x ◦ ∈ R n + , and let C x ◦ as the associated stoichiometric compatibility class. Assume there exists an RLF Lyapunovfunction and suppose that x ( t ) is bounded, Then the equilibrium set E x ◦ is Lyapunov stable. If, in addition, V satisfies the LaSalle’s Condition, then x ( t ) → E x ◦ as t → ∞ (i.e., thepoint to set distance of x ( t ) to E x ◦ tends to ). Furthermore, any isolated equilibriumrelative to C x ◦ is asymptotically stable. If V satisfies the LaSalle’s Condition, and all the trajectories are bounded, then: if thereexists x ∗ ∈ E x ◦ , which is isolated relative to C x ◦ then it is unique, i.e., E x ◦ = { x ∗ } .Furthermore, it is a globally asymptotically stable equilibrium relative to C x ◦ . Remark 5.
Note that the RLF considered can not be used to establish boundedness of solutions,as it may fail to be proper. Therefore, we need to resort to other methods so that boundednesscan be guaranteed. For instance, if the network is conservative, i.e the exists w ∈ R n + such that w T Γ = 0 , which ensures the compactness of C x ◦ .3.2. Uncertain Systems Framework: Reaction Coordinates In this subsection the function W R,x e is assigned to be R , as in the case of PWLR functions.As arbitrary monotone kinetics are allowed in our formulation of the CRN family N A,B , thesystem (5) with kinetics K A can be viewed as an uncertain system. However, this system is notin the form of the traditional types of parameter uncertainties treated in the literature. In thissubsection, we show that shifting the analysis of the system to reaction coordinates enables toview it as a linear parameter varying (LPV) system for which existence of a common Lyapunovfunction directly yields a robust Lyapunov function for the original CRN.Let r ( t ) := R ( x ( t )) , then we have: ˙ r ( t ) = ∂R∂x ( x ( t ))Γ r ( t ) = ρ ( t )Γ r ( t ) , (8)where ρ ( t ) := ∂R∂x ( x ( t )) . We can write ρ ( t ) as a conic combination of individual partial deriva-tives as follows: ∂R∂x ( x ( t )) = ρ ( t ) = X i,j : α ij > ρ ji ( t ) E ji , (9)where [ ρ ( t )] ji = ρ ji ( t ) , and [ E ji ] j ′ i ′ = 1 if ( j ′ , i ′ ) = ( j, i ) and zero otherwise.Let s denote the number of elements in the support of ∂R/∂x , and let κ : { , .., s } → { ( i, j ) : α ij > } be an indexing map. Then, we can write (8) as: ˙ r = X i,j : α ij > ρ ji ( t ) E ji Γ r = s X ℓ =1 ρ ℓ ( t )Γ ℓ r, (10)where Γ ℓ = e j γ Ti , ρ ℓ ( t ) = ρ ji ( t ) , with ( i, j ) = κ ( ℓ ) , and { e j } νj =1 denotes the canonical basis of R ν . Hence, equation (10) represents a linear parameter-varying system which has s nonnegative7ime-varying parameters { ρ ( t ) , .., ρ s ( t ) } and the system matrix belongs to the conic hull of theset of rank-one matrices { Γ , ..., Γ s } .Hence, we have the following definition. Definition 3 (Common Lyapunov Function) . A function ˜ V : ¯ R ν + → ¯ R + is said to be a Lyapunovfunction for the linear system ˙ r = Γ ℓ r if it is locally Lipschitz, nonnegative, has a negativesemi-definite time-derivative along the trajectories of the linear system, and ker ˜ V ⊂ ker Γ ℓ .Furthermore, ˜ V is said to be a common Lyapunov function for the set of linear systems { ˙ r =Γ r, ..., ˙ r = Γ s r } if it is a Lyapunov function for each of them, and ker ˜ V = T sℓ =1 ker Γ ℓ . Equipped with the above definitions, we are ready to state the main result for this Section. Itsproof is deferred to the appendix for the sake of readability.
Theorem 2 (Equivalence btw. CLF and RLF in reaction coordinates) . Given the system (5) .There exists a common Lyapunov function ˜ V : ¯ R ν + → ¯ R + for the set of linear systems { ˙ r =Γ r, ..., ˙ r = Γ s r } if and only if ( ˜ V , R ) induces the Robust Lyapunov function parameterized as V R ( x ) = ˜ V ( R ( x )) for the CRN family N A,B . Remark 6.
Since the zero matrix belongs to the conic hull of { Γ , ..., Γ s } , asymptotic stabil-ity can’t be established by the mere existence of the common Lyapunov function. A LaSalle’sargument is needed as will be mentioned in the following Section. Remark 7.
Although the dynamics of concentrations (5) or the dynamics of the extent of reaction (14) define positive systems, the differential linear inclusion that can be defined from (8) is nota positive linear differential inclusion. More precisely, this means that it is not necessarily truethat all the matrices that belong to conic hull { Γ , ..., Γ s } are Metzler.3.3. Dual Robust Lyapunov Function: Species Coordinates The RLF introduced in the previous subsection is a function of R ( x ) . We investigate nowRLFs that are functions of the difference x − x e . This is carried out in a manner that is dual towhat has been done in the previous subsection.In order to present the dual framework, an alternative representation of the system dynamicscan be adopted. Consider a CRN as in (5), and let x e be an equilibrium. Then, there exists x ′′ ( x ) ∈ ¯ R n + such that (5) can written equivalently as: ˙ x = Γ ∂R∂x ( x ′′ )( x − x e ) , x (0) ∈ C x e (11)The existence of x ′′ := x e + ε x ( x − x e ) for some ε x ∈ [0 , follows by applying theMean-Value Theorem to R ( x ) along the segment joining x e and x .8et z = x − x e , then similar to the previous section, the conic combination (9) can be usedto rewrite (11) as: ˙ z = Γ ∂R∂x ( x ′′ )( x − x e ) = s X ℓ =1 ρ ℓ ( t )Γ E ℓ z = s X ℓ =1 ρ ℓ ( t )Γ i ℓ e Tj ℓ z, (12)where ρ ℓ ( t ) = ∂R jℓ ∂x iℓ ( x ′′ ( x ( t )) , and Γ i is the i th column of Γ . Therefore, the system dynamicshas been embedded in the linear differential inclusion with vertices { Γ i e Tj , ..., Γ i s e Tj s } .Let D T be a matrix with columns that are the basis vectors of ker Γ T . The following theoremcan be stated. Theorem 3 (Equivalence CLF and RLF, species coordinates) . Given the system (5) . There existsa common Lyapunov function ˆ V : R n → ¯ R + for the set of linear systems { ˙ z = (Γ i e Tj ) z, ..., ˙ z =(Γ i s e Tj s ) z } , on the invariant subspace { z : D T z = 0 } if and only if ( ˜ V , W x e ) , W x e = x − x e induces the Robust Lyapunov function parameterized as V x e ( x ) = ˆ V ( x − x e ) for the CRN family N A,B .Proof.
Let V ( x ) = ˆ V ( x − x e ) , then when ∂ ˆ V /∂z exists we can write: ˙ V = ∂ ˆ V∂z ˙ z = ∂ ˆ V∂z s X ℓ =1 ρ ℓ ( t )Γ i ℓ e Tj ℓ z = s X ℓ =1 ρ ℓ ( t ) ∂ ˆ V∂z Γ i ℓ e Tj ℓ z ! . Since we have assumed that ˆ V is a common Lyapunov function for the set of linear systems { ˙ z = (Γ i e Tj ) z, ..., ˙ z = (Γ i s e Tj s ) z } the proof can proceed in both directions in a similar wayto the proof of Theorem 2. Notice that the constraint D T z = 0 is needed since D T ˙ x ( t ) ≡ isimplicit in the structure of the original system (5). We show next that if ˜ V used in reaction coordinates satisfies a relatively mild additionalassumption, then the Lyapunov function of the form ˆ V ( x − x e ) can be used, where x e is anequilibrium point for (5).To this end, the following theorem can be stated. Theorem 4 (Converting Lyapunov functions between reaction and species coordinates) . Let V ( x ) = ˜ V ( R ( x )) represent an RLF for the network family N A,B . If there exists ˆ V : R n → ¯ R + such that for all r ∈ R ν : ˜ V ( r ) = ˆ V (Γ r ) , (13) then V ( x ) = ˆ V ( x − x e ) represents an RLF for the network family N A,B , where x e is anyequilibrium point for (5) . roof. Condition 1 in Definition 1 is clearly satisfied. It remains to show the second condition.Let z = x − x e . Then, whenever ˆ V is differentiable: ˙ V ( x ) = ∂ ˆ V ( x − x e ) ∂z ˙ x = ∂ ˆ V ( x − x e ) ∂z Γ R ( x ) , Before proceeding, we prove two statements: First, from (13), we get ( ∂ ˜ V ( r ) /∂r ) = ( ∂ ˆ V (Γ r ) /∂z )Γ .Second, note that x − x e ∈ Im (Γ) , hence there exists R ∈ R ν such that Γ R = x − x e , where R can always be chosen nonnegative by assumption AS. Hence, where ˆ V is differentiable, we canuse (11) to write: ˙ V ( x ) = ∂ ˆ V ( x − x e ) ∂r Γ ∂R ( x ′′ ) ∂x ( x − x e ) = ∂ ˜ V ( R ) ∂r ∂R ( x ′′ ) ∂x Γ R = s X ℓ =1 ρ ℓ ∂ ˜ V ( R ) ∂r Γ ℓ R ≤ , where the last inequality follows from (29). Lemma 17 implies that ˙ V ( x ) ≤ for all x . Remark 8. V has a simpler structure than V since it depends on x − x e . However, it canbe noted in the proof that for each specific choice of x e , the Lyapunov function ˆ V ( x − x e ) isnonincreasing only along solutions contained in C x e . Relationship to the Extent of Reaction Formulation.
Recall that the extent of reaction [19] is defined as: ξ ( t ) = R t R ( x ( τ )) dτ + ξ (0) . If x ( t ) ∈ C x e ,then ∃ ξ ∗ ≥ such that x e − x ◦ = Γ ξ ∗ . We set ξ (0) := ξ ∗ . Hence, we can write: Γ ξ ( t ) = x ( t ) − x e , (14)and ˙ ξ = R ( x e + Γ ξ ) , ξ (0) := ξ ∗ , (15)which is the extent-of-reaction ODE representation of the dynamics of the CRN.Therefore, we state the following result. Corollary 5.
Let Γ be given and ˜ V : ¯ R ν + → ¯ R + be a function. Assume there exists ˆ V : R n + → ¯ R + such that for all r ∈ R ν , ˜ V ( r ) = ˆ V (Γ r ) . If ˜ V is a common Lyapunov functionfor { ˙ r = e j γ Ti r, ..., ˙ r = e j s γ Ti s r } where ( i ℓ , j ℓ ) = κ ( ℓ ) for all ℓ ∈ { , , . . . , s } , then ˜ V ( ξ ) isnonnegative, and nonincreasing along the trajectories of ˙ ξ = R ( x e + Γ ξ ) for any R ∈ K A .Proof. The first two statements follow from Theorems 2 and 4. We prove now the third state-ment. Using (13), (14) we get ˜ V ( ξ ) = ˆ V ( x − x e ) . Therefore, the required statement followsfrom the result that ˆ V ( x − x e ) is nonincreasing along the trajectories of (5).10 . Application to PWLR Lyapunov Functions In the previous papers [10, 11] the concept of Piecewise Linear in Rate (PWLR) Lyapunovfunctions has been introduced based on a direct analysis of the CRN. Such functions satisfy theconditions of Definition 1, and hence they are Robust Lyapunov functions. In this section weshow that those results can be interpreted in the uncertain systems framework introduced above.This also allows to provide alternative algorithms for the existence and construction of PWLRfunctions.Consider a CRN (5) with a Γ ∈ R n × r . Two representation of the PWLR Lyapunov functionhave been discussed. Given a partitioning matrix H ∈ R p × r such that ker H = ker Γ . PWLRLyapunov functions are piecewise linear in rates, i.e., they have the form: V ( x ) = ˜ V ( R ( x )) ,where ˜ V : R ν → R is a continuous PWL function given as ˜ V ( r ) = | c Tk r | , r ∈ ±W k , k = 1 , .., m/ , where the regions W k = { r ∈ R ν : Σ k Hr ≥ } , k = 1 , .., m form a proper conic partition of R ν , while { Σ k } mk =1 are signature matrices with the property Σ k = − Σ m +1 − k , k = 1 , .., m/ .The coefficient vectors of each linear component can be collected in a matrix C = [ c , .., c m ] T ∈ R m × r . If the function ˜ V is convex, then we have the following simplified representation of V : V ( x ) = k CR ( x ) k ∞ . This representation reminds of the ℓ ∞ -norm Lyapunov functions that have been used for linearsystems in [14]. In fact, the next theorem establishes the link between the results introduced in[11] for checking candidate PWLR functions based on direct analysis and previous work on ℓ ∞ Lyapunov functions using the framework introduced in the previous section.
Proposition 6.
Given Γ and H . Let V = ˜ V ◦ R be a candidate continuous nonnegative PWLRwith C = [ c ... c m ] T ∈ R m × r . Then ( ˜ V , R ) induces an RLF if and only if: ker C = ker Γ , and there exists matrices { Λ ℓ } sℓ =1 ⊂ R m × m such that Λ ℓ H = − C Γ ℓ , (16) and λ ℓk Σ k > , where Λ ℓ = [ λ ℓ T ...λ ℓm/ T ] T .If ˜ V is convex, then the second condition can be replaced with ) there exists Metzler matrices { Λ ℓ } sℓ =1 ⊂ R m × m such that Λ ℓ ˜ C = ˜ C Γ ℓ , (17) and Λ ℓ = 0 for all ℓ = 1 , .., s , where ˜ C = [ C T − C T ] T .Proof. The proof can be carried out by performing elementary algebraic manipulations on theresults presented in [11, Theorems 4,5]. The details are omitted for the sake of space.
Remark 9.
The symmetries in equation (17) imply that it can be written equivalently as: C Γ ℓ = ˜Λ ℓ C, (18) where ˜Λ ℓ is an m × m matrix which is defined by subtracting the upper m × m blocks of Λ ℓ from each other. ˜Λ ℓ satisfies: max k ˜ λ ( ℓ ) kk + X j = k | ˜ λ ( ℓ ) kj | ≤ . (19) This is exactly the condition that ℓ ∞ -norm Lyapunov functions need to satisfy for a linear system[13, 20]. This shows that Theorem 2 provides the framework to utilize the existing linear stabilityanalysis techniques in the literature to construct robust Lyapunov functions for nonlinear systemssuch as CRNs. For example, we can verify ℓ Lyapunov functions of the form V ( x ) = k CR ( x ) k directly by replacing condition (19) by max k ˜ λ ( ℓ ) kk + X j = k | ˜ λ ( ℓ ) jk | ≤ , (20) instead of converting them to the ℓ ∞ -norm form. Three construction algorithms have presented in [11] and we revisit the second one here.Before proceeding to the construction algorithm, we need to introduce the concept of a neighborto a region. Fix k ∈ { , .., m/ } . Consider H : for any pair of linearly dependent rows h Ti , h Ti eliminate h Ti . Denote the resulting matrix by ˜ H ∈ R ˜ p × ν , and let ˜Σ , .., ˜Σ m the correspondingsignature matrices. Therefore, the region can be represented as W k = { r | ˜Σ k ˜ Hr ≥ } . The distance d r between two regions W k , W j is defined to be the Hamming distance between ˜Σ k , ˜Σ j .Hence, the set of neighbors of a region W k are defined as: N k = { j ∈ { , , . . . , m } : d r ( W j , W k ) = 1 } , W k is one which differs only by the switching ofone inequality. Denote the index of the switched inequality by the map s k ( . ) : N k → { , .., p } .For simplicity, we use the notation s kℓ := s k ( ℓ ) .We use Theorem 2 to show that the problem of constructing a PWLR Lyapunov functionover a given partition, i.e. a given H , can be solved via linear programming. However, insteadof encoding the nondecreasingness condition into precomputed sign patterns as in the previouschapter, we use here alternative conditions which are stated in the following proposition. Proposition 7.
Given the system (5) and a partitioning matrix H ∈ R p × r . Consider the linearprogram: Find c k , ξ k , ζ k ∈ R ν , Λ ℓ ∈ R m × m , η kj ∈ R ,k = 1 , .., m ; j ∈ N k , ℓ = 1 , .., s, subject to c Tk = ξ Tk Σ k H,C Γ ℓ = − Λ ℓ H, λ ℓk Σ k ≥ ,c k − c j = η kj σ ks kj h s kj ,ξ k ≥ , T ξ k > , Λ ℓ ≥ . Then there exists a PWLR RLF with partitioning matrix H if and only if there exists a feasiblesolution to the above linear program that satisfies ker C = ker Γ satisfied. Furthermore, thePWLR RLF can be made convex by adding the constraints η kj ≥ . Remark 10.
A natural choice for H is H := Γ . The Lyapunov function reduces then to: V ( x ) = k diag ( ξ k ) ˙ x k , R ( x ) ∈ W k , which generalizes the Lyapunov function presented in [6]. Remark 11.
The LaSalle’s Condition can be verified via a graphical algorithm described in § III-F in [11].4.2. The Dual PWL Lyapunov Function In § III-C it has been shown that if there exists ˆ V such that ˜ V ( r ) = ˆ V (Γ r ) , then there exists adual RLF for the same network family. In the case of PWLR Lyapunov functions condition 1 inProposition 6 implies that this condition is always fulfilled. Hence, consider a PWLR Lyapunovfunction defined with a partitioning matrix H as in (4.1). By Proposition 6 and the assumptionthat ker H = ker Γ , there exists G ∈ R p × n and B ∈ R m × n such that H = G Γ and C = B Γ .Similar to {W} mk =1 , we can define the regions: V k = { z | Σ k Gz ≥ } , k = 1 , .., m, V k has nonempty interior iff W k has nonempty interior.Therefore, as the pair ( C, H ) specify the PWLR function fully, also the pair ( B, G ) specifies thefunction: ˆ V ( z ) = b Tk z, when Σ k Gz ≥ , where B = [ b , ..., b m ] T . If ˜ V is convex, then it can be written in the form: V ( x ) = k CR ( x ) k ∞ .Similarly, the convexity of ˆ V implies that V ( x ) = k B ( x − x e ) k ∞ , where the latter is the Lyapunov function used in [12].Theorem 4 established that if ˜ V ( R ( x )) is an RLF, then ˆ V ( x − x e ) is an RLF also. Thefollowing theorem shows that converse holds also for PWLR RLFs, however, it is worth recallingfrom Remark 8 that ˜ V ( R ( x )) is nonincreasing for all initials conditions, while ˆ V ( x − x e ) isnonincreasing only on C x e . Theorem 8.
Given (5) . Then, if there exists G ∈ R p × n and B ∈ R m × n such that: ( B Γ , G Γ) defines a PWLR RLF, then ( B, G ) defines a dual PWL RLF. ( B, G ) defines a dual PWL RLF, then ( B Γ , G Γ) defines a PWLR RLF. The proof is presented in the appendix.
Remark 12.
Since D T ( x − x e ) = 0 for x ∈ C x e , then if k B ( x − x e ) k ∞ is an RLF, then k ( B + Y D T )( x − x e ) k ∞ is also an RLF for an arbitrary matrix Y . Furthermore, since Theorem 8 hasshown that the reaction-based and the species-based representations are equivalent; it is easierto check and construct RLFs in the reaction-based formulation and they hold the advantage ofbeing decreasing over all stoichiometry classes.
5. Properties of Graphically Stable Networks
It has been shown in [11] that the Jacobian of any network admitting a PWLR RLF is P ,which implies that all principal minors are nonnegative. This can be used to show the followingresult. Theorem 9 (Robust Nonsingularity) . Given (5) . Assume that it is a graphically stable network.If for some realization R ∈ K A there exists a point in the interior of a proper stoichiometric classsuch that the reduced Jacobian is non-singular at it, then the reduced Jacobian is non-singularin the interior of R n + for any realization of the kinetics. This implies that any positive equilibriumof this network is isolated and non-degenerate relative to its class. roof. Recall that for a GSN the negative Jacobian is P for any choice of R ∈ K A [11]. Usingthe Cauchy-Binet formula [21], let I ⊂ { , .., n } be an arbitrary subset so that | I | = k . Thecorresponding principal minor can be written as: det I (cid:18) − Γ ∂R∂x (cid:19) = X J ⊂{ ,..,ν } , | J | = k det( − Γ IJ ) det (cid:18) ∂R∂x JI (cid:19) ( ∗ ) = X ι a ι Y ℓ ∈ L ι ⊂{ ,..,s } ρ ℓ , where ( ∗ ) refers to the fact that the sum can be expressed as a linear combination of productsof ρ , ..., ρ s . We claim that the coefficients a ι are all nonnegative. To show this, assume forthe sake of contradiction that there is some a ι ∗ negative. If we set all ρ ’s to zero except theones appearing in the ι th ∗ term, then this implies that the corresponding principal minor can benegative; a contradiction.Now, the theorem can be proven by noting that the reduced Jacobian is non-singular iff thesum of all k × k principal minors of the negative Jacobian is positive, where k = rank (Γ) . Sinceit is assumed that there exists a point for which the reduced Jacobian is non-singular, this impliesthat the sum of principal minors is positive for some choice of ρ , .., ρ s . Since all of the principalminors are nonnegative, then at least one of them is positive. By AK4, that principal minor stayspositive for any choice of positive ρ , .., ρ s , i.e. it stays positive over the interior of R n + . An important result links injectivity of a map to the notion of P -matrix, [22]. This states thata map is injective if its Jacobian is a P matrix. For an injective vector-field, if an equilibriumexists, then it is unique. This notion has been studied extensively for reaction networks, andspecialized by means of appropriate graphical conditions, [21, 23].Hence, the following theorem follows: Theorem 10 (Uniqueness of Positive Equilibria of Graphically Stable Network) . If N A,B isGS, then it can not admit multiple nondegenerate positive equilibria in a single stoichiometriccompatibility class. Furthermore, if there exists an isolated non-degenerate positive equilibrium x e , relative to C x e , then it is unique.Proof. It has been established in [23, Appendix B] that the network can not admit multiplenondegenerate positive equilibria in a single stoichiometric compatibility class if the Jacobian is P . Hence, the first statement follows.For the second statement, Theorem 9 has shown that the the existence of an non-degeneratepositive equilibrium x e ensures that the reduced Jacobian is non-singular on the interior of theorthant. In order to show uniqueness, assume for the sake of contradiction that there exists15 = x e , y ∈ C x e such that Γ R ( y ) = 0 . Then the fundamental theorem of calculus implies, R ( x e ) − Γ R ( y ) = Γ Z ∂R∂x ( tx e + (1 − t ) y ) ( x e − y ) dt = Γ ∂R∂x ( x ∗ )( x e − y ) , where x ∗ = t ∗ x e + (1 − t ∗ ) y , and t ∗ ∈ (0 , . The existence of t ∗ is implied by the integralmean-value theorem. Since x ∗ ∈ C ◦ x e , then the reduced Jacobian at x ∗ is non-singular relative toIm Γ . Since x e − y ∈ Im Γ , then y = x e : a contradiction. Remark 13.
Since the Jacobian is P , then if arbitrary inflows and outflow were added to ev-ery species of a GSN the resulting Jacobian would be a P matrix [24]; this is also known asthe continuous-flow stirred tank reactor (CFSTR) version of the network. The CFSTR networkis injective. This shows that our framework has a direct relationship to the recent results oninjectivity [25], P matrices for reaction networks [21, 23] and concordance [26].5.3. Persistence An important dynamical property in the context of positive systems is persistence . Informallyit holds if any solution initialized in the interior of the positive orthant will not asymptoticallyapproach its boundary. In [17] is shown that if a conservative network does not have criticalsiphons, then it is persistent. Note that this is a graphical property which is independent of thespecific realization of the kinetics involved.In this section it is shown that persistence can be established for conservative GSNs undersuitable conditions detailed in the following Theorem.
Theorem 11 (Absence of Types of Critical Siphons) . Given (5) . Consider the network family N A,B . Assume there exists a critical siphon P , and let Λ( P ) be the set of output reactions of P .Then, N A,B is a not GS, i.e., the network does not admit a PWLR Lyapunov function if any ofthe following conditions is satisfied. P is a critical deadlock. the network is conservative and for some realization of the network family there exists apoint in the interior of a proper stoichiometric compatibility class on which the reducedJacobian is nonsingular, the network is conservative and ker Γ is one-dimensional. The following theorem follows immediately from [17][Theorem 2] and Theorem 11.
Theorem 12 (Persistence of A Class of GSNs) . Given (5) . Assume the network is conservative.Then the network family N A,B is persistent if ker Γ is one-dimensional, or by removing the reverse of some reactions the reduced stoichiometry matrix has a one-dimensional kernel and the associated network is GS, or there exists an non-degenerate positive equilibrium in the interior of some proper stoichio-metric class for a realization of N A,B . The second item in Theorem 12 follows by the fact that the inclusion of a reverse of a reactiondoes not create a critical siphon.
We have shown that the existence of a PWLR Lyapunov function implies that it is a commonLyapunov function for all linear systems that belong to a linear differential inclusion.In fact, one of the properties of systems that admits piecewise linear Lyapunov function isthat a stable equilibrium can not have purely imaginary eigenvalues [27]. Hence, the reduced Ja-cobian at a non-degenerate equilibrium can not admit pure imaginary eigenvalues which impliesthe following Theorem:
Theorem 13 (Exponential Stability) . Given (5) that admits a PWLR function. If a positiveequilibrium x e is non-degenerate relative to C x e , then it is exponentially asymptotically stable. Remark 14.
For a conservative network the state space is compact. Therefore, the existence ofa non-degenerate equilibrium implies that it is unique (Theorem 10), it is locally exponantiallystable (Theorem 13), the level sets of the Lyapunov function are always invariant (Theorem 1)and the network is persistent (Theorem 12). However, global asymptotic stability (GAS) can notbe claimed directly without a LaSalle argument. It is potentially possible that the equilibrium isnot GAS; for instance, there can be a limit cycle living in the boundary of the basin of attractionthat attracts the outside trajectories. Despite the fact that this seems unlikely, it can not be pre-cluded completely without a proof. Therefore, the graphical algorithm for verifying the LaSalle’sargument presented in [11] is still needed to claim GAS.
6. Relationship to Contraction Analysis
Contraction analysis is an approach to stability investigation focused on the relative behaviourof solutions, rather than on their deviations from a nominal trajectory (such as an equilibriumpoint). This area of research is as old as the concept of contraction mappping, however, it hassparked growing interest in the control systems community in relationship to the analysis ofdynamical systems [28], [29].The are several formulations of contraction theory. We are going to present the formulationthat utilizes matrix measures, or logarithmic norms.17 efinition 4 (Logarithmic Norms) . For a given induced matrix norm k . k ∗ on R n × n , the associ-ated matrix measure (or logarithmic norm) can be defined as follows for a matrix A ∈ R n × n : µ ∗ ( A ) := lim sup h → + k I + hA k ∗ − h . (21) Note that the same definition applies if k . k ∗ is a semi-norm. Remark 15.
The logarithmic norm can be evaluated for the standard norms. For instance, thefollowing expression can be used for the ℓ ∞ norm: µ ∞ ( A ) = max i a ii + X j = i | a ij | . (22) Note that this expression is identical to the one appearing in (19) . For a dynamical system, negativity of the logarithmic norm can be linked to contraction. Thisresult has been stated in different forms, refer to the tutorial [29] for more details. We state theresult as follows.
Theorem 14 ([29]) . Consider a dynamical system ˙ x = f ( x ) defined on a convex subset X of R n . Let | · | ∗ be a norm in R n and k . k ∗ the induced matrix norm on R n × n . Assume that ∀ x ∈ X, µ ∗ (cid:18) ∂f∂x ( x ) (cid:19) ≤ c. Then for any two solutions ϕ ( t ; x ) , ϕ ( t ; y ) of the dynamical system, the following conditionholds: | ϕ ( t ; x ) − ϕ ( t ; y ) | ∗ ≤ e ct | ϕ (0; x ) − ϕ (0; y ) | ∗ . (23)Note that if c < the solutions of the system are exponentially contracting. If c = 0 , thenthe system is non-expansive. The choice of the norm plays a crucial role even with respect to di-agonal weighings. The result above have been applied to CRNs before by weighting the ℓ -normwith a diagonal matrix [30]. The link between the logarithmic norms and norm-based Lyapunovfunctions has been established before [20]. Since convex PWLR Lyapunov functions are ( ∞ ) -norms weighted by a non-square matrix, a similar result can be expected to hold. We state in thefollowing Theorem the precise relationship between convex PWLR functions introduced beforeand contraction analysis. Theorem 15 (Relationship to Contraction Analysis) . Given the extent of reaction representation of CRNs ˙ ξ = R ( x e + Γ ξ ) . Assume that thereexists a convex PWLR function ˜ V ( ξ ) = k Cξ k ∞ , and let µ C be the logarithmic norm ssociated. Let the associated Jacobian be: J ( ξ ) := ∂R∂x Γ . Then, ∀ ξ, µ C ( J ( ξ )) ≤ . Hence, the system is non-expansive on the subspace (ker Γ) ⊥ , where R ν = ker Γ ⊕ (ker Γ) ⊥ . Given the ODE ˙ x = Γ R ( x ) . Assume that there exists convex PWL function V ( x ) = k B ( x − x e ) k ∞ , and let µ B be the logarithmic norm associated. Let the associated Jaco-bian be: J ( x ) := Γ ∂R∂x . Then, ∀ x ∈ C x e , µ B ( J ( x )) ≤ . Hence, the system is non-expansive in each stoichiometric class C x e . Remark 16.
The upper bound µ ∞ (cid:16)P sℓ =1 ρ ℓ ˜Λ ℓ (cid:17) can be identical to zero especially if m ≥ n .Therefore, an analogous concept to a LaSalle argument need to be introduced. This is discussedin the next section.6.1. Variational Dynamics and LaSalle Argument for Contraction Analysis In a recent work, Forni and Sepulchre [31] have proposed a Lyapunov framework for contrac-tion analysis using a so-called
Finsler structure . In order to minimize the background needed,we apply it directly to our context. The Finsler-Lyapunov function for the system (15) is V F : T ¯ R ν + → ¯ R + . If contraction analysis was carried out with respect to the norm: k . k ∗ : ξ Cξ k ∞ , then the corresponding Finsler-Lyapunov function would be V F ( δξ ) = k Cδξ k ∞ . The Finsler structure is given by the mapping δξ
7→ k
Cδξ k ∞ . Therefore, V F can be consideredas a Lyapunov function for the variational system: ˙ δξ = ∂R∂x Γ δ ξ. (24)The Finsler structure induces a distance function, which is, in this case, d F ( x, y ) = k C ( x − y ) k ∞ . Hence, if V F is strictly decreasing then this would imply that this system is incrementallyasymptotically stable with respect to the distance function [31], which is equivalent to the resultgiven by Theorem 14. However, if strict decreasingness does not hold, the Finsler-Lyapunovframework for contraction analysis has the advantage of accommodating a LaSalle’s invarianceprinciple that can be used to show strict contraction. In fact, the same algorithm proposed in [11]19an be used for the Finsler-Lyapunov function to test the LaSalle’s Condition. The followingtheorem states the result. Theorem 16 (Strict Contraction) . Given (5) . Assume that V ( x ) = k CR ( x ) k ∞ is a PWLRLyapunov function that satisfies Proposition 7 in [11] . Then the trajectories of (15) are expo-nentially contractive with respect to the norm k . k ∗ : ξ
7→ k Cξ k ∞ in directions orthogonal to ker Γ .Proof. As per [31, Theorem 2], we need to show that if a trajectory lives in ker ˙ V F , then it is anequilibrium for the variational system (24). Since V F is a piecewise function it can be studiedper partition regions as before. Hence let V F ( δξ ) = c Tk δξ when δξ ∈ W k . Then, ˙ V F ( δξ ) = c Tk ˙ δξ = c Tk ∂R∂x Γ δξ. Note that the expression above is analogous to (25) in [11], where sgn (Γ δξ ) can be made constantin each partition region. Therefore, the arguments of the proof of [11, Proposition 7] can bereplicated to show that the algorithm proposed in [11] implies that when δξ ( t ) ∈ ker ˙ V F for all t ≥ , then Γ δξ ( t ) ≡ . Remark 17.
Parallel to the duality expounded in the previous sections, the results can also beinterpreted for the variational system in species coordinates: ˙ δx = Γ ∂R ( x ) ∂x δx, Dδx (0) = 0 . (25) Then, ˆ V ( δx ) is a Lyapunov function for (25) .
7. Conclusions
We have presented a theoretical complement for [11]. It has been shown that PWLR Lya-punov functions introduced in [11] can be interpreted as defining a common Lyapunov functionfor a linear differential inclusion in the reaction coordinates. Furthermore, a dual Lyapunovfunction in the species coordinates has been defined. Many additional properties of graphicallystable networks have been shown. Examples of biochemical networks for which our frameworkis applicable can be found in [32].
Appendix: Proofs
Before proving the results of the paper, we need to state and prove the following Lemma:20 emma 17.
Let ˙ x := f ( x ) , and let V : ¯ R n + → ¯ R + be a locally Lipschitz function such that: ∂V ( x ) ∂x f ( x ) ≤ , whenever ∂V ( x ) ∂x exists , then ˙ V ( x ) ≤ for all x .Proof. Since V is assumed to be locally Lipschitz, Rademacher’s Theorem implies that it isdifferentiable (i.e., gradient exists) almost everywhere [33]. Recall that for a locally Lipschitzfunction the Clarke’s gradient at x can be written as ∂ C V ( x ) := co ∂V ( x ) , where: ∂V ( x ) := n p ∈ R n : ∃ x i → x with ∂V ( x i ) /∂x exists, such that , p = lim i →∞ ∂V ( x i ) /∂x o . Let p ∈ ∂V ( x ) and let { x i } ∞ i =1 be the corresponding sequence. By the assumption stated in thelemma, ( ∂V ( x i ) /∂x ) f ( x i ) ≤ , for all i . Hence, the definition of p implies that p T f ( x ) ≤ .Since p was arbitrary, the inequality holds for all p ∈ ∂V ( x ) .Now, let p ∈ ¯ ∂V ( x ) where p = P i λ i p i is a convex combination of any p , ..., p n +1 ∈ ∂V ( x ) .By the inequality above, p T f ( x ) = P i λ i ( p Ti f ( x )) ≤ . Hence, p T f ( x ) ≤ for all p ∈ ¯ ∂V ( x ) .As in [33], the Clarke’s derivative of V at x in the direction of f ( x ) can be written as D Cf ( x ) V ( x ) =max { p T f ( x ) : p ∈ ¯ ∂V ( x ) } . By the above inequality, we get D Cf ( x ) V ( x ) ≤ for all x . Sincethe Dini’s derivative is upper bounded by the Clarke’s derivative, we finally get: ˙ V ( x ) := lim sup h → + V ( x + hf ( x )) − V ( x ) h ≤ lim sup h → + y → x V ( y + hf ( x )) − V ( y ) h =: D Cf ( x ) V ( x ) ≤ , for all x . Proof of Theorem 2.
We show the existence of the common Lyapunov function implies theexistence of the RLF. Nonnegativity of V follows from the nonnegativity of ˜ V . Let ( i, j ) = κ ( ℓ ) ,recall that Γ ℓ = e j γ Ti , hence ker ˜ V = T sℓ =1 ker Γ ℓ = ker Γ . Therefore, R ( x ) ∈ ker V iff Γ R ( x ) = 0 , which establishes the positive-definiteness of V .We assumed that ˜ V has a negative semi-definite time-derivative for every linear system in theconsidered set. Hence, when ˜ V is differentiable, we can write ( ∂ ˜ V /∂r )Γ ℓ r ≤ , ℓ = 1 , ..., s .Hence, for any ρ , ..., ρ s ∈ ¯ R + : s X ℓ =1 ρ ℓ ∂ ˜ V∂r Γ ℓ r ≤ , when ( ∂V ( r ) /∂r ) exists . (26)21herefore, when ˜ V is differentiable: ˙ V ( x ) = ∂ ˜ V∂R ∂R∂x ( x )Γ R ( x ) = ∂ ˜ V∂R X i,j : α ij > ∂R j ∂x i ( x ) E ji Γ R ( x ) , (27)where ∂ ˜ V /∂R := ( ∂ ˜ V /∂r ) (cid:12)(cid:12)(cid:12) r = R ( x ) .Now, denote ρ ℓ = ∂R j ∂x i ( x ) , which is nonnegative by A3. This allows us to write: ˙ V ( x ) = s X ℓ =1 ρ ℓ ∂ ˜ V∂R E ji Γ R ( x ) (28) = s X ℓ =1 ρ ℓ ∂ ˜ V∂R Γ ℓ R ( x ) ≤ , for almost all x. (29)The last inequality follows from (26). Using Lemma 17, ˙ V ( x ) ≤ for all x , and for all R ∈ K A .In order to show the other direction, almost all properties outlined in Definition 3 are clearlysatisfied, we just show nonincreasingness. Assume that there exists ℓ such that ˜ V ( r ) is notnonincreasing along the trajectories of ˙ r = Γ ℓ r . Consider the corresponding term in (28). Since V ( R ( x )) is a Lyapunov function for any choice of admissible rate reaction function R , choose ρ ℓ = ∂R j ∂x i to be large enough such that ˙ V ( x ) ≥ for some x ; a contradiction. (cid:4) Proof of Theorem 8.
The first statement follows from Theorem 4. In order to show the secondstatement, let V ( x ) = b Tk ( x − x e ) , for x − x e ∈ V k . We will show that V ( x ) = c Tk R ( x ) , for R ( x ) ∈ W k is nondecreasing. Without loss of generality, the partition matrix can be written inthe form: G = [ I ˆ G T ] T . This representation implies that the sign of x − x e is determined inevery region V k . Now, assume that x − x e ∈ V ◦ k , then: ˙ V ( x ) = b Tk Γ R ( x ) = c Tk R ( x ) ≤ c Tk R ( x e ) , for all R ∈ K A . Let R j ( x ) ∈ supp R ( x ) , and let α ij > . Since R is nondecreasing by A3, if sgn ( x i − x e i ) sgn ( c kj ) > , there exists R ∈ K A such that ˙ V ( x ) ≥ . Hence, this implies that theinequality sgn ( c kj ) sgn ( x i − x e i ) ≤ holds. Fix j , if there exists i , i such that α i j , α i j > and sgn ( x i − x e i ) sgn ( x i − x e i ) < , then σ kj := 0 . Otherwise, σ kj := sgn ( x i − x e i ) forsome i such that α ij > .Hence, in order to have ˙ V ( x ) ≤ for all R ∈ K A we need that σ kj ( x i − x e i ) ≥ whenever x − x e ∈ V k , for all k, j, i with α ij > . By Farkas Lemma [34], this is equivelant to theexistence of λ kji ∈ ¯ R p + , ζ kji ∈ R ι : σ kj e Ti = λ Tkji Σ k G + ζ Tkji D, (30)22here D T ∈ R ι × n is a matrix whose columns are basis vectors for ker Γ T .If we multiply both sides of (30) by Γ from the left, then we get condition C4 in [11, Theorem 4]which necessary and sufficient for ˙ V ( x ) = ddt ( c Tk R ( x )) ≤ . (cid:4) Proof of Theorem 11.
Assume P is a critical siphon for the petri-net associated with Γ , and let n p = | P | . Let Λ( P ) be the set of output reactions of P , and let ν p = | Λ( P ) | .Before we prove item 1 of Theorem 11, the following lemma is needed. Lemma 18.
Consider a network family N A,B . Let P be a set of species that does not containthe support of a conservation law; let its indices be numbered as { , ..., n p } . Then, there existsa nonempty-interior region { r | Σ k Γ r ≥ } with a signature matrix Σ k that satisfies σ k = ... = σ kn p = 1 .Proof. Assume the contrary. This implies that ∩ n p i =1 { R | γ Ti R > } T ∩ ni = n p +1 { R | σ i γ Ti R > } = ∅ for all possible choices of signs σ i = ± . However, R r can be partitioned into a unionof all possible half-spaces of the form ∩ ni = n p +1 { R | σ i γ Ti R ≥ } . Therefore, this implies that ∩ n p i =1 { R | γ Ti R > } = ∅ . By Farkas Lemma, this implies that there exists λ ∈ R t satisfying λ > such that [ λ T . Therefore, P contains the support of the conservation law [ λ T T ;a contradiction.Therefore, we can state the proof of the first item: Proof of Theorem 11-1).
Without loss of generality, let { , ..., n p } be the indices of the speciesin P . Using Lemma 18, there exists a nonempty-interior sign region S k , ≤ k ≤ m s with asignature matrix Σ k that satisfies σ k = ... = σ kn p = 1 . Since Λ( P ) = R , this implies b kj ≤ for all j = 1 , .., ν . However, this is not allowable by [11, Theorem 9] since ζ Tk B k v ≤ for all v ∈ ker Γ ∩ ¯ R n + and for any choice of admissible ζ k .In order to proceed, an existence result of equilibria is needed: Lemma 19.
Consider a network family N A,B . Let P be a critical siphon, and let Ψ P be theassociated face. If the network is conservative, then for any proper stoichiometric compatibility C , there exists an equilibrium x e of (5) such that x e ∈ Ψ P ∩ C .Proof. The set Ψ P ∩ C is compact, forward invariant, and convex, since the both sets Ψ P , C are as such. Hence, the statement of the lemma follows directly from the application of theBrouwer’s fixed point theorem on the associated flow.We are ready now to prove the second item of Theorem 11. Proof of Theorem 11-2).
By Lemma 19, there exists an equilibrium in Ψ P . Since it is assumedthat there exists an isolated equilibrium in interior, Theorem 10 implies that N A,B is not GS.23efore concluding the proof, a simple lemma is stated and proved:
Lemma 20.
Let x e be an equilibrium of (5) . Let ˜ P be a set of species that correspond to { , .., n }\ supp ( x e ) . Then, ˜ P is a siphon.Proof. Assume that ˜ P is not a siphon, then there exists some X i ∈ ˜ P and R j ∈ R such that X i is a product of R j and R j = Λ( ˜ P ) . At the given equilibrium, all negative terms in the expressionof ˙ x i vanish since x ei = 0 . Since X i is not a reactant in R j this implies β ij > , α ij = 0 then R j ( x ) has a strictly positive coefficient which implies ˙ x i > ; a contradiction.Hence, we are ready to conclude the proof of Theorem 11: Proof of Theorem 11-3).
By Lemma 19, there exists an equilibrium x ∗ ∈ Ψ P such that Γ R ( x ∗ ) =0 . Since dim(ker Γ) = 1 , this implies that R ( x ∗ ) = tv for some t ≥ . Consider the case t = 0 .This implies R ( x ∗ ) = 0 . Then, P ⊂ ˜ P := { , .., n }\ supp ( x ∗ ) . ˜ P is a siphon by Lemma 20,and since P ⊂ ˜ P it is a critical deadlock. However, by Theorem 11-1), N A,B is not GS. If t > ,this implies that P = ∅ ; a contradiction. Proof of Theorem 14.
Write the logarithmic norm expression using (21): µ C ( J ( ξ )) = lim sup h → + h (cid:18)(cid:13)(cid:13)(cid:13)(cid:13) I + h ∂R∂x Γ (cid:13)(cid:13)(cid:13)(cid:13) C − (cid:19) The expression above includes the induced matrix norm. Using the definition of the inducedmatrix norm we proceed as follows: (cid:13)(cid:13)(cid:13)(cid:13) I + h ∂R∂x Γ ξ (cid:13)(cid:13)(cid:13)(cid:13) C = sup k Cξ k ∞ =1 (cid:13)(cid:13)(cid:13)(cid:13) C (cid:18) I + h ∂R∂x Γ (cid:19) ξ (cid:13)(cid:13)(cid:13)(cid:13) ∞ ( ⋆ ) = sup k Cξ k ∞ =1 (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) Cξ + h s X ℓ =1 ρ ℓ ( x ) Ce j γ Ti ξ (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∞ ( ♣ ) = sup k Cξ k ∞ =1 (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) Cξ + h s X ℓ =1 ρ ℓ ˜Λ ℓ Cξ (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∞ ≤ sup k Cξ k ∞ =1 (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) I + h s X ℓ =1 ρ ℓ ˜Λ ℓ (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∞ k Cξ k ∞ (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) I + h s X ℓ =1 ρ ℓ ˜Λ ℓ (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∞ , where ( ⋆ ) is by (9) and ( ♣ ) is by (18). Therefore, the expression of the logarithmic norm abovecan be written as: µ C ( J ( ξ )) ≤ µ ∞ s X ℓ =1 ρ ℓ ˜Λ ℓ ! ≤ s X ℓ =1 ρ ℓ µ ∞ (˜Λ ℓ ) = 0 , (31)where the inequalities follow by the subadditivity of the logarithmic norm and (19).Note since C has nonempty kernel space, then k Cξ k ∞ is a semi-norm. However, Theorem14 requires a norm. This can be remedied by studying the system in directions orthogonal to24 er Γ by defining a transformation of coordinates using a matrix T = [ ˆ T , v , .., v ν − ν r ] T , where { v , .., v ν − ν r } is a basis of ker Γ , ν r = rank (Γ) , and ˆ T is chosen so that T is invertible.Defining ˆ ξ = T ξ , then the first ν r coordinates of ˆ ξ are decoupled from the rest. Hence, inequality(31) can be established similarly for the reduced Jacobian which is the upper right ν r × ν r blockof T J ( ξ ) T − . The norm in the reduced subspace is k . k ˆ C : z
7→ k ˆ Cz k ∞ , z ∈ R ν , where ˆ C is m × ν r defined as the nonzero columns of CT − . The equations Ce j γ Ti = ˜Λ ℓ C are equivalent to ( CT − )( T e j γ Ti T − ) = ˜Λ ℓ ( CT − ) . Therefore, everything goes through in the reduced space,and the same upper bound is valid.The same argument can be replicated for to prove the second item in the theorem. Thisis accomplished by utilizing the alternative representation (11), the rank-one decomposition ofthe Jacobian (12), and noting that conditions (18) can be written as: B Γ i e Tj = ˜Λ ℓ B + Y ℓ D T for some matrices Y , .., Y ℓ . The reduced space argument can be carried out also by using atransformation matrix T = [ ˆ T , D ] T . (cid:4) ReferencesReferences [1] E. D. Sontag, Structure and stability of certain chemical networks and applications to the kinetic proofreadingmodel of T-cell receptor signal transduction, IEEE Transactions on Automatic Control 46 (7) (2001) 1028–1047.[2] D. Angeli, A tutorial on chemical reaction network dynamics, European Journal of Control 15 (3-4) (2009) 398–406.[3] J. E. Bailey, Complex biology with no parameters, Nature Biotechnology 19 (6) (2001) 503–504.[4] F. Horn, R. Jackson, General mass action kinetics, Archive for Rational Mechanics and Analysis 47 (2) (1972)81–116.[5] M. Feinberg, The existence and uniqueness of steady states for a class of chemical reaction networks, Archive forRational Mechanics and Analysis 132 (4) (1995) 311–370.[6] H. Maeda, S. Kodama, Y. Ohta, Asymptotic behavior of nonlinear compartmental systems: nonoscillation andstability, IEEE Transactions on Circuits and Systems 25 (6) (1978) 372–378.[7] D. Angeli, P. De Leenheer, E. Sontag, Graph-theoretic characterizations of monotonicity of chemical networks inreaction coordinates, Journal of mathematical biology 61 (4) (2010) 581–616.[8] F. Blanchini, E. Franco, Structurally robust biological networks, BMC systems biology 5 (1) (2011) 74.[9] M. Ali Al-Radhawi, D. Angeli, Lyapunov functions for the stability of a class of chemical reaction networks, in: the 20th International Symposium on Mathematical Theory of Networks and Systems , Melbourne, Australia, 2012.URL: .[10] M. Ali Al-Radhawi, D. Angeli, Piecewise linear in rates Lyapunov functions for complex reaction networks, in:Proceedings of the 52nd IEEE Control and Decision Conference (CDC), 2013, pp. 4595–4600.[11] M. A. Al-Radhawi, D. Angeli, New approach for the stability of complex reaction networks: Piecewise linear inrates Lyapunov functions, IEEE Transactions on Automatic Control 61 (1) (2016) 76–89.[12] F. Blanchini, G. Giordano, Piecewise-linear Lyapunov functions for structural stability of biochemical networks,Automatica 50 (10) (2014) 2482 – 2493.[13] A. P. Molchanov, E. S. Pyatnitskii, Lyapunov functions that specify necessary and sufficient conditions of absolutestability of nonlinear nonstationary control systems. I,III, Automation and Remote Control 47 (1986) 344–354,620–630.
14] A. Polanski, On infinity norms as Lyapunov functions for linear systems, IEEE Transactions on Automatic Control40 (7) (1995) 1270–1274.[15] F. Blanchini, Nonquadratic Lyapunov functions for robust control, Automatica 31 (3) (1995) 451–461.[16] M. Ali Al-Radhawi, D. Angeli, Robust Lyapunov functions for complex reaction networks: An uncertain systemframework, in: Proceedings of the IEEE 53rd Conference on Decision and Control (CDC), 2014, pp. 3101–3106.[17] D. Angeli, P. De Leenheer, E. D. Sontag, A Petri net approach to persistence analysis in chemical reaction networks,in: I. Queinnec, S. Tarbouriech, G. Garcia, S. Niculescu (Eds.), Biology and Control Theory: Current Challenges,Springer, 2007, pp. 181–216.[18] T. Yoshizawa, Stability theory by Liapunov’s second method, Mathematical Society of Japan, Tokyo, 1966.[19] B. L. Clarke, Stability of complex reaction networks, in: I. Prigogine, S. Rice (Eds.), Advances in ChemicalPhysics, Volume 43, John Wiley & Sons, 1980, pp. 1–215.[20] H. Kiendl, J. Adamy, P. Stelzner, Vector norms as Lyapunov functions for linear systems, IEEE Transactions onAutomatic Control 37 (6) (1992) 839–842.[21] M. Banaji, P. Donnell, S. Baigent, P matrix properties, injectivity, and stability in chemical reaction systems, SIAMJournal on Applied Mathematics 67 (6) (2007) 1523–1547.[22] D. Gale, H. Nikaido, The jacobian matrix and global univalence of mappings, Mathematische Annalen 159 (2)(1965) 81–93.[23] M. Banaji, G. Craciun, Graph-theoretic approaches to injectivity and multiple equilibria in systems of interactingelements, Communications in Mathematical Sciences 7 (4) (2009) 867–900.[24] A. Berman, R. J. Plemmons, Nonnegative matrices in the mathematical sciences, Academic Press, New York, 1979.[25] G. Craciun, M. Feinberg, Multiple equilibria in complex chemical reaction networks: I. The injectivity property,SIAM Journal on Applied Mathematics (2005) 1526–1546.[26] G. Shinar, M. Feinberg, Concordant chemical reaction networks, Mathematical biosciences 240 (2012) 92–113.[27] E. B. Castelan, J. C. Hennet, On invariant polyhedra of continuous-time linear systems, in: Proceedings of the 30thIEEE Conference on Decision and Control, 1991, pp. 1736–1741.[28] W. Lohmiller, J.-J. Slotine, On contraction analysis for non-linear systems, Automatica 34 (6) (1998) 683–696.[29] Z. Aminzare, E. D. Sontag, Contraction methods for nonlinear systems: A brief introduction and some open prob-lems, in: Proc. 53rd IEEE Conf. on Decision and Control, 2014, pp. 3835–3847.[30] G. Russo, M. Di Bernardo, E. D. Sontag, Global entrainment of transcriptional systems to periodic inputs, PLoScomputational biology 6 (4) (2010) e1000739.[31] F. Forni, R. Sepulchre, A differential Lyapunov framework for contraction analysis, IEEE Transactions on Auto-matic Control 59 (3) (2014) 614–628.[32] M. Ali Al-Radhawi, D. Angeli, Construction of robust Lyapunov functions for reaction networks, in: Proceedingsof the 15th European Control Conference (ECC), 2016, pp. 928–935.[33] F. H. Clarke, Y. Ledyaev, R. Stern, P. Wolenski, Nonsmooth analysis and control theory, Springer, New York, 1997.[34] R. T. Rockafellar, Convex analysis, Princeton University Press, New Jersey, 1970.14] A. Polanski, On infinity norms as Lyapunov functions for linear systems, IEEE Transactions on Automatic Control40 (7) (1995) 1270–1274.[15] F. Blanchini, Nonquadratic Lyapunov functions for robust control, Automatica 31 (3) (1995) 451–461.[16] M. Ali Al-Radhawi, D. Angeli, Robust Lyapunov functions for complex reaction networks: An uncertain systemframework, in: Proceedings of the IEEE 53rd Conference on Decision and Control (CDC), 2014, pp. 3101–3106.[17] D. Angeli, P. De Leenheer, E. D. Sontag, A Petri net approach to persistence analysis in chemical reaction networks,in: I. Queinnec, S. Tarbouriech, G. Garcia, S. Niculescu (Eds.), Biology and Control Theory: Current Challenges,Springer, 2007, pp. 181–216.[18] T. Yoshizawa, Stability theory by Liapunov’s second method, Mathematical Society of Japan, Tokyo, 1966.[19] B. L. Clarke, Stability of complex reaction networks, in: I. Prigogine, S. Rice (Eds.), Advances in ChemicalPhysics, Volume 43, John Wiley & Sons, 1980, pp. 1–215.[20] H. Kiendl, J. Adamy, P. Stelzner, Vector norms as Lyapunov functions for linear systems, IEEE Transactions onAutomatic Control 37 (6) (1992) 839–842.[21] M. Banaji, P. Donnell, S. Baigent, P matrix properties, injectivity, and stability in chemical reaction systems, SIAMJournal on Applied Mathematics 67 (6) (2007) 1523–1547.[22] D. Gale, H. Nikaido, The jacobian matrix and global univalence of mappings, Mathematische Annalen 159 (2)(1965) 81–93.[23] M. Banaji, G. Craciun, Graph-theoretic approaches to injectivity and multiple equilibria in systems of interactingelements, Communications in Mathematical Sciences 7 (4) (2009) 867–900.[24] A. Berman, R. J. Plemmons, Nonnegative matrices in the mathematical sciences, Academic Press, New York, 1979.[25] G. Craciun, M. Feinberg, Multiple equilibria in complex chemical reaction networks: I. The injectivity property,SIAM Journal on Applied Mathematics (2005) 1526–1546.[26] G. Shinar, M. Feinberg, Concordant chemical reaction networks, Mathematical biosciences 240 (2012) 92–113.[27] E. B. Castelan, J. C. Hennet, On invariant polyhedra of continuous-time linear systems, in: Proceedings of the 30thIEEE Conference on Decision and Control, 1991, pp. 1736–1741.[28] W. Lohmiller, J.-J. Slotine, On contraction analysis for non-linear systems, Automatica 34 (6) (1998) 683–696.[29] Z. Aminzare, E. D. Sontag, Contraction methods for nonlinear systems: A brief introduction and some open prob-lems, in: Proc. 53rd IEEE Conf. on Decision and Control, 2014, pp. 3835–3847.[30] G. Russo, M. Di Bernardo, E. D. Sontag, Global entrainment of transcriptional systems to periodic inputs, PLoScomputational biology 6 (4) (2010) e1000739.[31] F. Forni, R. Sepulchre, A differential Lyapunov framework for contraction analysis, IEEE Transactions on Auto-matic Control 59 (3) (2014) 614–628.[32] M. Ali Al-Radhawi, D. Angeli, Construction of robust Lyapunov functions for reaction networks, in: Proceedingsof the 15th European Control Conference (ECC), 2016, pp. 928–935.[33] F. H. Clarke, Y. Ledyaev, R. Stern, P. Wolenski, Nonsmooth analysis and control theory, Springer, New York, 1997.[34] R. T. Rockafellar, Convex analysis, Princeton University Press, New Jersey, 1970.