Chemical Reaction Network Decomposition Technique for Stability Analysis
CChemical Reaction Network Decomposition Technique for Stability Analysis ∗ Yafei Lu † , Chuanhou Gao † , and Denis Dochain ‡ Abstract.
This paper develops the concept of decomposition for chemical reaction networks, based on which anetwork decomposition technique is proposed to capture the stability of large-scale networks charac-terized by a high number of species, high dimension, high deficiency, and/or non-weakly reversiblestructure. We present some sufficient conditions to capture the stability of a network (may possessany dimension, any deficiency, and/or any topological structure) when it can be decomposed intoa complex balanced subnetwork and a few 1-dimensional subnetworks (and/or a few two-speciessubnetworks), especially in the case when there are shared species in different subnetworks. Theresults cover encouraging applications on autocatalytic networks with some frequently-encounteredbiochemical reactions examples of interest, such as the autophosphorylation of PAK1 and Aurora Bkinase, autocatalytic cycles originating from metabolism, etc.
Key words. chemical reaction networks, stability, autocatalytic, mass-action system, decomposition
AMS subject classifications.
1. Introduction.
Chemical reaction networks (CRNs) have a pivotal role in various fieldsincluding chemistry, system biology, molecular biology as well as other potential applications.The research on CRNs aims at studying the connection between the networks’ structure andunderlying dynamical characteristics, e.g. oscillations, persistence, stability, etc. Thereinto,the issue of stability for CRNs, especially in the mass-action kinetics setting, has receivedconsiderable critical attention, which is also the concern in the current paper.The study refers to the stability for CRNs with general structure is difficult in generalwhile a number of researchers [3, 6, 9–12, 27] have reported that certain networks equippedwith special structure (e.g. reversible, or weakly reversible) exhibit convergence behaviors.Horn and Jackson [18] proposed the well-known deficiency zero theorem, which stated that aweakly reversible deficiency zero mass action system (MAS) is complex balanced and has aunique equilibrium in every positive stoichiometric compatibility class. Meanwhile, each suchequilibrium is proved to be asymptotically stable by taking the pseudo-Helmholtz functionas a Lyapunov function. Later, this result was extended by Feinberg [11], named as thedeficiency one theorem, in which the deficiency can be relaxed, i.e., not necessarily equal tozero. Another way to reach the stability of equilibria for complex balanced MASs (also containsdetailed balanced MASs) was introduced in [28], which established a compact mathematicalformulation to exhibit the dynamics of the network through the complex graph. Further, thestudies of [19, 20] have applied the results on complex balanced MASs to more CRNs. It says ∗ Submitted to the editors. .
Funding:
This work was funded by the National Nature Science Foundation of China under Grant No. 12071428,11671418, and the Zhejiang Provincial Natural Science Foundation of China under Grant No. LZ20A010002. † School of Mathematical Sciences, Zhejiang University ([email protected] (Y. Lu), [email protected] (C.Gao, Correspondence)). ‡ ICTEAM, UCLouvain, Bˆatiment Euler, avenue Georges Lemaˆıtre 4-6, 1348 Louvain-la-Neuve, Belgium ([email protected] (D. Dochain). This manuscript is for review purposes only. a r X i v : . [ m a t h . D S ] F e b Y. LU, C. GAO AND D. DOCHAIN that if an MAS can be transformed into a complex balanced MAS through the linear conjugacyapproach, then the former shares the same stability as the latter. Several other attempts havebeen made to the stability for more CRNs by constructing their corresponding Lyapunovfunctions. For instance, the scaling limits of non-equilibrium potential [2] were taken as theLyapunov function for birth-death processes from a microscopic perspective; the piecewiselinear in rates Lyapunov function [1] was put forward for some balanced MASs; the paper [22]developed a generalized pseudo-Helmholtz function to the stability analysis for a generalbalanced MAS by taking advantage of the reconstruction and reverse reconstruction strategies.In the meanwhile, the generalized pseudo-Helmholtz function also succeeded in settling thestability for complex-balanced-produced MASs [30]. Recently, Fang and Gao [8] proposed asystematic approach called Lyapunov function PDEs to the stability for general CRNs, whichhas been validated well on complex balanced, 1-dimensional, special high dimensional andcomplex balanced produced CRNs so far.However, except for the above special cases, in most cases, such as the biological systems,many CRNs are non-weakly reversible, of any deficiency, and of any arbitrary dimension.More importantly, when the dimensions of the networks become large, the associated stabilityproblem becomes challenging. Especially, as the rapid development of systems biology andsynthetic biology, studying the dynamical behaviors of complex and large biochemical reac-tions has drawn great attention. Some investigators focuses on the properties of large-scaleCRNs by joining small networks (e.g., cross-talk) or decomposing networks (e.g., inhibition),like identifiability [13], multistationarity, stability [25] and stationary distribution [15]. Theymainly use the information of the small parts to infer the properties of the large networksynthesized by these parts. In particular, these contributions will promote the research oncross-talk [7, 14, 24], which refers to the situation where two or more signaling pathways withthe same components affect each other, and shows great application potential in medicine.In line with the above mentioned studies, this paper is devoted to deriving the stability oflarge-scale CRNs by decomposing them into low dimensional networks and/or parts havingspecial structures (e.g. reversible, weakly reversible). Note that the networks we are concernedwith are all non-weakly reversible, of any dimension, and of any deficiency. To be specific,we first show the results that if the decomposition of a CRN contains a complex balancedCRN and several 1-dimensional/two-species CRNs, then we can determine its stability byutilizing the information given by subnetworks. i.e. the corresponding Lyapunov functionsfor the subnetworks. These results cover two basic cases (i) the subnetworks have disjointspecies sets, and (ii) the subnetworks have shared species. Further, such a dimension reduc-tion decomposition approach is applied to analyze the stability for some practical biologicalreactions, including the autophosphorylation of PAK1 and Aurora B kinase, autocatalytic cy-cles, etc. Moreover, we deal with a common class of autocatalytic reactions as an application.It proves that the stability for such CRNs can be captured in terms of the characteristicsof the small parts: two-species autocatalytic CRNs generated from decomposition. In somesense, the present study may offer some significant insights into the field of stability analysisto synthetic biology and artificial life.The remaining part of the paper proceeds as follows. Section 2 states the motivation ofthis study. In section 3, we discuss the definition of decomposition and the decompositiontechnique that serves for capturing the stability of a CRN. Section 4 suggests a kind of special
This manuscript is for review purposes only.
RN DECOMPOSITION FOR STABILITY 3 decomposition that includes two-species subnetworks, and the sufficient condition to renderstability is also presented in this special case. In section 5, we consider some applications ofthe proposed methods to autocatalytic CRNs. Finally, section 6 concludes the paper.
Notation: R n , R n ≥ , R n> : n -dimensional real space, non-negative real space, positive real space, respec-tively. Z n ≥ : n -dimensional non-negative integer space. x v · i : x v · i = (cid:81) dj =1 x v ji j , where x ∈ R d , v · i ∈ Z d and 0 = 1.Ln( x ) : Ln( x ) = (ln x , · · · , ln x n ) (cid:62) , where x ∈ R n> . C i ( · ; ∗ ) : the set of i th continuous differentiable functions from ” · ” to ”*”. s.t. : such that. e i : unit vector with i th element being one while others being zero.
2. Motivation statement.
A CRN N = ( S , C , R ) equipped with mass-action kinetics iscalled an MAS, denoted by M = ( S , C , R , K ), where S is the species set, C the complex set, R the reaction set, and K the reaction rate constant set. More terminologies about CRNsand MASs can be found in Appendix 1: CRNs and related stability resutlts. CRNs or MASscan be grouped according to the dimension of stoichoimetric subspace into 1-dimensional,2-dimensional networks, etc., or according to structure into weakly reversible and non-weaklyreversible networks. For an M = ( S , C , R , K ) with the dynamics of(1) d x d t = ΓΞ( x ) = k i x v .i , x ∈ R n ≥ , the stability issue in the Lyapunov sense is about how to construct an available Lyapunovfunction suggesting the equilibrium x ∗ constrained by ΓΞ( x ∗ ) = 0 stable. It was reported thatthere have been available Lyapunov functions for complex balanced networks [18] (i.e., theclassical pseudo-Helmholtz free energy function (A.5)) and 1-dimensional networks [8] (i.e.,(A.6)).Although it becomes relatively routine to generate Lyapunov functions for the above twoclasses of networks, there are no systematic methods for others to construct the Lyapunovfunctions, especially for those extremely complicated ones (or large-scale) characterized byhigh number of species, high dimension, high deficiency, and/or non-special structure, likemetabolic regulatory networks, gene regulatory networks in biological systems. It will pose agreat challenge to perform the stability analysis for them. However if we observe those large-scale networks from the viewpoint of part or locality, the subnetwork may be very special, evento be weakly reversible and/or 1-dimensional. We use the following example [15] to illustrateour idea. This manuscript is for review purposes only.
Y. LU, C. GAO AND D. DOCHAIN
Example 1.
The network route takes S (cid:47) S (cid:111) (cid:47) S (cid:111) (cid:47) S (cid:111) ,S + S (cid:47) (cid:47) S , S + S (cid:47) (cid:47) S , (2) 2 S (cid:47) S + S (cid:111) , S (cid:47) S (cid:111) , S (cid:47) (cid:47) S (cid:122) (cid:122) S + S (cid:79) (cid:79) , which looks rather complicated, and has dimension , deficiency and non-weakly reversiblestructure. However, its subnetworks, as an example of a kind of decomposition, may includea weakly reversible network on { S , S , S } given by N : 3 S (cid:47) S (cid:111) , S (cid:47) (cid:47) S (cid:122) (cid:122) S + S (cid:79) (cid:79) , and three -dimensional CRNs on { S , S } , { S , S } , { S , S } successively, N : S (cid:47) S (cid:111) , S + S (cid:47) (cid:47) S , N : S (cid:47) S (cid:111) , S + S (cid:47) (cid:47) S , N : S (cid:47) S (cid:111) , S (cid:47) S + S (cid:111) . If every subnetwork is stable (relatively easy to be checked), is it possible for the originalnetwork to be stable too, or can the Lyapunov function for the original network be derivedfrom those known ones for special subnetworks? This naive idea motivates us to analyze thestability of a large-scale or complicated network through decomposing it into some small-scaleand/or simple subnetworks for which the Lyapunov functions are known.
3. Decomposition technique.
There are two possibilities in decomposing a network intosome subnetworks: one is that all subnetworks have no disjoint species, the other is that thereare shared species among subnetworks. We discuss these two cases in this section, but givethe definition of decomposition first.
Definition 1. (Decomposition).
For an M = ( S , C , R , K ) given by (1) that admits anequilibrium x ∗ ∈ R n> , if there are finitely many M ( p ) = ( S ( p ) , C ( p ) , R ( p ) , K ( p ) ) satisfying (1) S = ∪ S ( p ) , C = ∪ C ( p ) , R = ∪ R ( p ) , K = ∪ K ( p ) ; (2) every M ( p ) = ( S ( p ) , C ( p ) , R ( p ) , K ( p ) ) admits an equilibrium x ∗ ( p ) ∈ R n p > , in which eachentry is the same as that in x ∗ if they represent the concentration of the same species in S ,then N = ( S , C , R ) is said to be decomposable. Based on the decomposition definition, we have the following property about the equilib-rium.
This manuscript is for review purposes only.
RN DECOMPOSITION FOR STABILITY 5
Property 1.
Consider an M = ( S , C , R , K ) governed by (1) , and x ∗ ∈ R n> is an equilibriumin M . Assume M can be decomposed into a complex balanced M (0) = ( S (0) , C (0) , R (0) , K (0) ) and (cid:96) M ( p ) = ( S ( p ) , C ( p ) , R ( p ) , K ( p ) ) , p = 1 , ..., (cid:96) . If ∀ p ∈ { , · · · , (cid:96) } , x ∗ ( p ) ∈ R n p > , as defined in Definition , is a reaction vector balanced equilibrium in M ( p ) , then x ∗ isa generalized balanced equilibrium in M .Proof. Let each entry of x ∗ (0) ∈ R n > equals to the entry of x ∗ pointing to the same speciesin S , then x ∗ (0) is a complex balanced equilibrium in M (0) . Based on Definition A.8 andDefinition 1, it is straightforward to reach the conclusion.Note that the network decomposition will not generate or lose species, complex and reac-tion. Since every subnetwork is relatively simple and smaller-scale compared with the originalnetwork, it is possible to analyze the stability of every subnetwork and accordingly infer thatof the original network. There are two cases when implementing network composition, one is ∀ p (cid:54) = q, S ( p ) ∩ S ( q ) = ∅ , the other is S ( p ) ∩ S ( q ) (cid:54) = ∅ . In the following, we will discuss these twosituations respectively, and moreover, the subnetworks are limited to be complex balancedand 1-dimensional for convenience of stability analysis. This case is relativelytrivial, especially when the subnetworks are complex balanced and 1-dimensional. In thefollowing we directly give a conclusion.
Theorem 2.
For an M = ( S , C , R , K ) governed by (1) , let x ∗ ∈ R n> be an equilibrium in M . Suppose that M can be decomposed into a complex balanced M (0) = ( S (0) , C (0) , R (0) , K (0) ) and (cid:96) M ( p ) = ( S ( p ) , C ( p ) , R ( p ) , K ( p ) ) , p = 1 , ..., (cid:96) , and moreover, ∀ p (cid:54) = q ∈{ , , ..., (cid:96) } , S ( p ) ∩S ( q ) = ∅ . If for every 1-dimensional subnetwork there is ω (cid:62) p ∂∂x ( p ) h ( x ∗ ( p ) , < , then M is locally asymptotically stable at x ∗ , where h ( x ∗ ( p ) , and ω p follow the definitionof (A.7) , and x ∗ ( p ) ∈ R n p > is an equilibrium in M ( p ) , p = 1 , ..., (cid:96), n + (cid:80) (cid:96)p =1 n p = n. Proof.
The result is similar to Theorem 27 in [8] where a composition network containinga complex balanced MAS and some 1-dimensional MASs is proved asymptotically stable, andwe thus omit it here.
The situation will become in-tractable if the subnetworks share species in network decomposition. In this case, the dynamicbehaviors of the subsystems will affect each other through the shared species. However, wemanage to derive a sufficient condition to suggest the asymptotic stability of a class of net-works that can be decomposed into a complex balanced subnetwork and a few 1-dimensionalsubnetworks but with species shared among subnetworks.
Theorem 3.
Given an M = ( S , C , R , K ) with the dynamics (1) and an equilibrium x ∗ ∈ R n> , assume it can be decomposed into a complex balanced M (0) and (cid:96) M ( p ) ’s, p = 1 , · · · , (cid:96) , and moreover, ∀ p, q ∈ { , · · · , (cid:96) } , E p := S (0) (cid:84) S ( p ) (cid:54) = ∅ , S ( p ) (cid:84) S ( q ) ∈ (cid:83) (cid:96)p =1 E p or ∅ , and every M ( p ) is reaction vector balanced. Then M is locally asymptotically stable at x ∗ if for every p ∈ { , ..., (cid:96) } the following conditions are true (1) for all S j ∈ E p in v ( p ) · l −→ v (cid:48) ( p ) · l there is v ( p ) · m −→ v (cid:48) ( p ) · m such that v ( p ) jm = v (cid:48) ( p ) jl , v (cid:48) ( p ) jm = v ( p ) jl ,and v (cid:48) ( p ) jm − v ( p ) jm = v ( p ) jl − v (cid:48) ( p ) jl = 1 ;This manuscript is for review purposes only. Y. LU, C. GAO AND D. DOCHAIN (2) denote ˜ x ( p ) = ( (cid:78) S j ∈S ( p ) \ E p x ( p ) j ) (cid:62) , L ω p = { l : v (cid:48) ( p ) · l − v ( p ) · l = ( (cid:78) S i ∈ E p , ω (cid:62) p ) (cid:62) } , and R ω p = { l : v (cid:48) ( p ) · l − v ( p ) · l = − ( (cid:78) S i ∈ E p , ω (cid:62) p ) (cid:62) } , it holds that ω (cid:62) p ∇ ˜ u p (˜ x ( p ) ) | ˜ x ( p ) =˜ x ∗ ( p ) > , (3) where ˜ u p (˜ x ( p ) ) = (cid:81) S i ∈ E p x ∗ i ( p ) (cid:80) R ωp k ( p ) l (cid:81) S j ∈S ( p ) \ E p x ( p ) j v ( p ) jl (cid:80) L ωp k ( p ) l (cid:81) S j ∈S ( p ) \ E p x ( p ) j v ( p ) jl . Proof.
The detailed proof can be found in Appendix 2. B. The proof of Theorem 3.
4. Special decomposition: a complex balanced subnetwork and some two-species sub-networks with shared species.
In this section, we weaken the conditions presented in The-orem 3 by considering the decomposition including two-species subnetworks, which are 1-dimensional essentially, instead of general 1-dimensional networks. Further, we extend theprevious decomposition method to deal with more complex CRNs for capturing stability.
Consider a class of two-species CRNs with dimen-sion 1 defined on the species set { S i , S j } , given by L ω : v · l k l (cid:47) (cid:47) v (cid:48)· l , R ω : v · l k l (cid:47) (cid:47) v (cid:48)· l , (4)where L ω = { l : v (cid:48)· l − v · l = ω (cid:62) = ( w i , w j ) } , and R ω = { l : v (cid:48)· l − v · l = − ω (cid:62) = − ( w i , w j ) } . Thisclass of CRNs satisfies that each v il = a ∀ l in L ω and each v jl = b ∀ l in R ω .The stability of the two-species MASs can be naturally reached by adopting the Lyapunovfunction given in (A.6) for all 1-dimensional MASs. However, we propose another one belowfor the same purpose, which yet will play an important role on weakening the conditions ofTheorem 3 when implementing the network decomposition technique. Lemma 4.
For a class of two-species CRNs defined above, suppose it is reaction vectorbalanced with an equilibrium x ∗ = ( x ∗ i , x ∗ j ) ∈ R > , then the function f ( x ) = − (cid:90) x i x ∗ i w − i ln t a c ij (cid:80) R ω k l t v il dt + (cid:90) x j x ∗ j w − j ln t b c ij (cid:80) L ω k l t v jl dt (5) with c ij = x ∗ ai (cid:80) Rω k l x ∗ vili = x ∗ bi (cid:80) Lω k l x ∗ vjlj can act as a Lyapunov function to suggest that x ∗ is locallyasymptotically stable if there are w − i (cid:88) R ω k l ( a − v il ) x ∗ vil − i < and w − j (cid:88) L ω k l ( b − v jl ) x ∗ vjl − j > . (7) This manuscript is for review purposes only.
RN DECOMPOSITION FOR STABILITY 7
Proof.
Firstly, we verify that f ( x ) is positive definite. The Hessian matrix of f ( x ) iscalculated by ∇ f ( x ) = diag (cid:32) − w − i (cid:80) R ω k l ( a − v il ) x v il − i (cid:80) R ω k l x v il i , w − j (cid:80) L ω k l ( b − v jl ) x v jl − j (cid:80) L ω k l x v jl j (cid:33) . Based on the continuity of ∇ f ( x ) with respect to x and conditions (6),(7), there exists aneighborhood of x ∗ , denoted by D ( x ∗ ), such that ∀ x ∈ D ( x ∗ ) (cid:84) S + ( x ∗ ), there is ∇ f ( x ) >
0, which means f ( x ) is locally strictly convex. This together with the fact ∇ f ( x ∗ ) = (cid:18) − w − i ln x ∗ ai c ij (cid:80) Rω k l x ∗ vili , w − j ln x ∗ bj c ij (cid:80) Lω k l x ∗ vjlj (cid:19) (cid:62) = (0 , (cid:62) yields f ( x ∗ ) = 0 and f ( x ) > x ∈ D ( x ∗ ) (cid:84) S + ( x ∗ ) but x (cid:54) = x ∗ .Next we prove that f ( x ) fulfills dissipativeness. The dynamics of the two-species MAS aregiven by (cid:40) ˙ x i = w i (cid:80) L ω k l x ai x v jl j − w i (cid:80) R ω k l x bj x v il i , ˙ x j = w j (cid:80) L ω k l x ai x v jl j − w j (cid:80) R ω k l x bj x v il i . (8)Then we have ˙ f ( x ) = ∇ (cid:62) f ( x ) ˙ x = − ln (cid:80) L ω k l x ai x v jl j (cid:80) R ω k l x bj x v il i (cid:18) (cid:88) L ω k l x ai x v jl j − (cid:88) R ω k l x bj x v il i (cid:19) ≤ . (9)Here, the last equality holds if and only if (cid:88) L ω k l x ai x v jl j − (cid:88) R ω k l x bj x v il i = 0 , which means x = x ∗ . Therefore, f ( x ) can be used as a Lyapunov function to capture the localasymptotic stability of x ∗ .In the following, we consider a special class of two-species networks, i.e., two-speciesautocatalytic CRNs. Autocatalytic CRNs is an active area of research on biological systems(e.g., DNA replications [26], population evolution, etc), life and artificial life [17, 21, 29] (todescribe self-organizing systems), and chemistry. The notion of autocatalytic reactions isusually used to describe a class of reactions that are catalyzed by the products. Definition 5 (a reduced version of autocatalytic CRNs [15]).
An MAS M = ( S , C , R , K ) issaid to be autocatalytic if it satisfies (1) all reactions have a net consumption of one S i and a net production of one S j , i.e., inthe form of S i + ( α j − S j k i,αj (cid:47) (cid:47) α j S j , α j ≥ , This manuscript is for review purposes only.
Y. LU, C. GAO AND D. DOCHAIN and for convenience, denote the collection of the above reactions by R i,j := { v · l −→ v (cid:48)· l ∈ R : v (cid:48)· l − v · l = e j − e i } , where e i and e j are unit vectors; (2) there must exist a pair of monomolecular reversible reactions; (3) if S i −→ S j ∈ R i,j ⊂ R and S l −→ S j ∈ R l,j ⊂ R , then the reactions in R i,j and R l,j contain reactions of the same molecularity such that there exists some c ∈ R > with c · ( k i, , ..., k i,n j ) = ( k l, , ..., k l,n j ) , where n j represents the highest integer of the reaction S i + ( n j − S j −→ n j S j . Any two-species autocatalytic CRN is a 1-dimensional network, and is thus locally asymp-totically stable at its equilibrium based on Lemma 4.
Corollary 6.
For any two-species autocatalytic M = ( S , C , R , K ) governed by Eq. (8) with L ω = R i,j , R ω = R j,i , w i = − , w j = 1 , assume it admits a reaction vector balanced equilib-rium x ∗ ∈ R > , the function f ( x ) = (cid:90) x i x ∗ i ln x i c ij (cid:80) R j,i k j,α i x α i − i dx i + (cid:90) x j x ∗ j ln x j c ij (cid:80) R i,j k i,α j x α j − j dx j (10) with c ij = x ∗ i (cid:80) R j,i k j,αi x ∗ αi − i = x ∗ j (cid:80) R i,j k i,αj x ∗ αj − j can behave as a Lyapunov function to show x ∗ is locally asymptotically stable if (cid:88) v · l −→ v (cid:48)· l ∈R i,j k i,α j (2 − α j ) x ∗ αj − j > and (cid:88) v · l −→ v (cid:48)· l ∈R j,i k j,α i (2 − α i ) x ∗ αi − i > are ture.Proof. By applying Lemma 4 to the two-species autocatalytic MAS, the result is obvious.
Example 2.
Given a two-species autocatalytic network with deficiency , which has thereaction route S k , (cid:47) S k , (cid:111) , (13) 2 S + S k , (cid:47) (cid:47) S ,S + 2 S k , (cid:47) (cid:47) S , S S + S k , (cid:111) (cid:111) k , (cid:47) (cid:47) S . This manuscript is for review purposes only.
RN DECOMPOSITION FOR STABILITY 9
The type of autocatalytic CRN can be used to model the foraging system in ants for the researchon decision-making problems [23]. Here, we focus on the properties of its equilibria. Bychoosing k , = k , = k , = 1 , k , = 3 , k , = 2 , k , = 4 , the dynamical equations are asfollows (cid:26) ˙ x = − (4 x + x x + x x ) + 2 x + x x + 3 x x , ˙ x = 4 x + x x + x x − (2 x + x x + 3 x x ) . We can calculate an equilibrium (1 , (cid:62) within the positive compatibility class { x + x = 2 } .According to Corollary , we can find a candidate Lyapunov function like Eq. (10) for thesystem f ( x ) = (cid:90) x ln 6 tt + 3 t + 2 dt + (cid:90) x ln 6 tt + t + 4 dt. Then it is easy to validate that (1 , (cid:62) satisfies Eqs. (6) and (7) , that is − x ∗ | x ∗ =1 = 1 > , and − x ∗ | x ∗ =1 = 3 > . Therefore, (1 , (cid:62) is locally asymptotically stable.Remark
7. If a two-species autocatalytic CRN is at-most-biomolecular (the sum of sto-ichiometric coefficients of any complex in the network is at most 2 [13]), i.e., every α i ≤ α j ≤
2, then this system must be locally asymptotically stable since Eqs. (6) and (7) aresatisfied naturally.
When the decomposition contains two-species networks instead of general 1-dimensional networks, we have the following result.
Theorem 8.
For an M = ( S , C , R , K ) described by (1) with an equilibrium x ∗ ∈ R n> ,assume it can be broken into a complex balanced M (0) = ( S (0) , C (0) , R (0) , K (0) ) and (cid:96) reac-tion vector balanced two-species M ( p ) = ( S ( p ) , C ( p ) , R ( p ) , K ( p ) ) , p = 1 , ..., (cid:96) , and moreover, ∀ p, S ( p ) = { S i , S j } ⊆ S and E p := S (0) (cid:84) S ( p ) (cid:54) = ∅ . Then the local asymptotic stability of x ∗ can be achieved if the following conditions are true, (1) let E = (cid:83) (cid:96)p =1 E p , for every S i ∈ E , the reaction in R w p in M ( p ) satisfies v ( p ) il − a ( p ) = w ip , (2) for any species S j ∈ S ( p ) (cid:84) S ( q ) , p, q ∈ { , · · · , (cid:96) } , the reactions in M ( p ) and M ( q ) contain the same stoichiometric coefficient of S j such that there is a constant c ∈ R > with ( k ( p )1 , · · · , k ( p ) r p ) = c · ( k ( q )1 , · · · , k ( q ) r q ) , r p = r q . (3) ∀ S j ∈ S\S (0) and ∀ p , there holds w − jp (cid:88) L ωp k ( p ) l ( b ( p ) − v ( p ) jl ) x ∗ v ( p ) jl − j > . (14) This manuscript is for review purposes only.
Proof.
The detailed proof can be found in Appendix 2. C. The proof of Theorem 8.
Remark
9. Unlike Theorem 3, the above result suggests that the decomposed two-speciesCRNs may share more species but not just the species included in S (0) .By combing Theorem 8 with Theorem 3, we can easily obtain another decompositionpattern that may be applied to more CRNs. Corollary 10.
Given a M = ( S , C , R , K ) governed by (1) with an equilibrium x ∗ ∈ R n> , as-sume it can be decomposed into a complex balanced M (0) and (cid:96) M ( p ) ’s satisfyingreaction vector balancing. Then the system is locally asymptotically stable at x ∗ if (1) ∀ p, q ∈ { , · · · , (cid:96) } , when M ( p ) , M ( q ) are both two-species CRNs and all their reactionssatisfy condition (2) in Theorem ; when M ( p ) , M ( q ) are two-species CRNs but do not satisfycondition (2) in Theorem , or they are not two species CRNs, there is S ( p ) (cid:84) S ( q ) ∈ (cid:83) (cid:96)p =1 E p or ∅ with E p := S (0) (cid:84) S ( p ) (cid:54) = ∅ ; (2) ∀ p = 1 , · · · , (cid:96) , for those M ( p ) s that are two-species CRNs and satisfy condition (2) inTheorem , conditions (1) and (3) in Theorem is satisfied, and for others conditions (1) and (2) in Theorem are satisfied. We revisit the network (2) in Example 1 to illustrate Theorem 8 and Corollary 10.
Example 3.
Assume this network is constrained by x + x + x = 3 , x + x = 2 and x + x = 2 . Under the rate constants given below it has an equilibrium x ∗ = (1 , , , , (cid:62) .For convenience of applying Corollary , we consider four subnetworks (1) a complex balanced M (0) on { S , S , S } with an equilibrium ( x ∗ , x ∗ , x ∗ ) (cid:62) = (1 , , (cid:62) ,i.e., N : 3 S (cid:47) S (cid:111) , S (cid:47) (cid:47) S (cid:122) (cid:122) S + S (cid:79) (cid:79) ;(2) two-species autocatalytic M (1) and M (2) on { S , S } and { S , S } , which own respec-tive reaction vector equilibrium ( x ∗ , x ∗ ) (cid:62) = (1 , (cid:62) and ( x ∗ , x ∗ ) (cid:62) = (1 , (cid:62) , N : S (cid:47) S (cid:111) , S + S (cid:47) (cid:47) S , and N : S (cid:47) S (cid:111) , S + S (cid:47) (cid:47) S ;(3) one -dimensional two-species M (3) on { S , S } , which has a reaction vector balancedequilibrium ( x ∗ , x ∗ ) (cid:62) = (1 , (cid:62) , N : S (cid:47) S (cid:111) , S (cid:47) S + S (cid:111) . This manuscript is for review purposes only.
RN DECOMPOSITION FOR STABILITY 11
From Definition and Property , the network is decomposable, and x ∗ = (1 , , , , (cid:62) is a generalized equilibrium in M . Further, we have E = { S } , E = { S } , E = { S } , S (1) (cid:84) S (3) = ∅ , S (2) (cid:84) S (3) = { S } ∈ E and S (1) (cid:84) S (2) = { S } / ∈ (cid:83) i =1 E i . Finally, we checkthe other needed conditions in Corollary . • condition (1) in Corollary holds since all the reactions in N and N satisfy condition (2) in Theorem . Although N is a two-species CRN, it does not meet condition (2) inTheorem . • condition (2) in Corollary holds: condition (1) in Theorem is true in M (1) ,since ω = ( − , , for S −→ S in R ω , v (1)12 − a (1) = 0 − ω , and in M (2) similaranalysis can be made to prove true; condition (3) in Theorem holds since M (1) and M (2) are two-species autocatalytic MASs and fulfill Eq. (14) . Besides, condition (1) in Theorem is true since in M (3) , S ∈ E , for S −→ S , S + S −→ S , there are S −→ S and S −→ S + S ; condition (2) in Theorem holds in M (3) : L ω = { l : v (cid:48) (3) · l − v (3) · l = (1 , − } and R ω p = { l : v (cid:48) (3) · l − v (3) · l = ( − , } , so it is not difficult to compute that ˜ u ( x ) = 2 + 2 x x + x , and w ∇ ˜ u | x ∗ =1 = 2 x ∗ + 4 x ∗ + 6(3 x ∗ + x ∗ ) = 34 > , which means that M (3) gratifies Eq. (3) , too.As a result, the equilibrium x ∗ of M is locally asymptotically stable, and the Lyapunov functionis given by f ( x ) = (cid:88) i =1 , , (1 − x i + x i ln x i ) + (cid:90) x ln 2 tt + 1 dt + (cid:90) x ln 3 t + t t dt. (15)
5. Applications: autocatalytic CRNs decomposed for stability.
In this section, the de-composition strategy is applied to autocatalytic networks for stability analysis.As stated in subsection 4.1, autocatalytic CRNs arise abundantly in the biochemical sys-tems, which is deemed as the fundamental to the life process (e.g, self-replications of RNAmolecules [16], metabolic pathways). A typical example is an autocatalytic metabolism, suchas the Calvin-Benson-Bassham cycle and glyoxylate cycle in central carbon metabolism [4].Thus the study of stability analysis for autocatalytic CRNs is helpful to understand and en-gineer the metabolic networks. The previous results in section 4 provide a feasible solution toaddress the issue of stability of autocatalytic CRNs.
Property 2.
Assume an autocatalytic M = ( S , C , R , K ) admitting an equilibrium x ∗ ∈ R n> can be decomposed into several two-species M ( p ) ’s ( p = 1 , · · · , (cid:96) ) according to every pair ( R i,j , R j,i ) . Then x ∗ is reaction vector balanced in M if and only if each x ∗ ( p ) = ( x ∗ i , x ∗ j ) ∈ R > , S i , S j ∈ S is reaction vector balanced in M ( p ) .Proof. The detailed proof can be found in Appendix 2. D. The proof of Property 2.
This manuscript is for review purposes only.
We thus can present a sufficient condition to suggest the stability of an autocatalyticnetwork through the decomposition strategy.
Theorem 11.
Given an autocatalytic M = ( S , C , R , K ) with an equilibrium x ∗ ∈ R n> ,suppose it can be decomposed into a few two-species autocatalytics M ( p ) ’s ( p = 1 , · · · , (cid:96) ) as-sociated with ( R i,j , R j,i ) . Then x ∗ is locally asymptotically stable if every M ( p ) is reactionvector balanced and gratifies Eqs. (6) and (7) .Proof. Combining Corollary 6 and Property 2, the result comes immediately.We give two examples to exhibit the validity of Theorem 11 and Property 2: one is aconstructed network, the other is a reduced metabolic autocatalytic network [4].
Example 4.
The following autocatalytic network is -dimensional and with deficiency , S (cid:47) S (cid:111) (cid:47) S (cid:111) (cid:47) S (cid:111) ,S + 2 S (cid:47) (cid:47) S , S + S (cid:47) (cid:47) S , (16) S + 2 S (cid:47) (cid:47) S , S + S (cid:47) (cid:47) S ,S + 2 S (cid:47) (cid:47) S , S + S (cid:47) (cid:47) S . Assume it obeys three conservation laws x + x = √ , x + x = √ and x + x = √ ,then it is easy to get that x ∗ = ( √ − , , , √ − ) (cid:62) is a reaction vector balanced equilibrium.We decompose the network into three two-species subsystems M ( p ) | p =1 with S (1) = { S , S } , S (2) = { S , S } , S (3) = { S , S } . It is also easy to verify that ( x ∗ , x ∗ ) (cid:62) = ( √ − , (cid:62) , ( x ∗ , x ∗ ) (cid:62) = ( √ − , (cid:62) and ( x ∗ , x ∗ ) (cid:62) = (1 , √ − ) (cid:62) are reaction vector balanced in M (1) , M (2) and M (3) , respectively. These results reveal that the MAS share the same equilibriumwith that of the subsystems, which is consistent with P roperty . Finally we can verify thatevery subsystem M ( p ) supports Eqs. (6) and (7) , so based on Theorem the MAS is locallyasymptotically stable at ( √ − , , , √ − ) (cid:62) . Example 5.
Consider an n -species autcatalytic MAS with dim S = n − and δ = n , S i k i, (cid:47) S i +1 k i +1 , (cid:111) , (17) S i + S i +1 k i, (cid:47) (cid:47) S i +1 ,i = 1 , · · · , n, S n +1 = S . The lower part in this network can be viewed as an autocatalytic cycle model [4] without inflowand outflow, which usually appears in metabolic networks.By selecting k i, = k i, = 1 , k i +1 , = 2 , i = 1 , · · · , n , and supposing this network to followthe mass conservation laws x i + x i +1 = 2 , i = 1 , · · · , n − , we get x ∗ = (1 , · · · , (cid:62) is a reactionThis manuscript is for review purposes only. RN DECOMPOSITION FOR STABILITY 13 vector balanced equilibrium in this MAS. Clearly, this network can be decomposed into n two-species autocatalytic CRNs according to reaction vectors, i.e., S (1) = ( S , S ) , S (2) = ( S , S ) , · · · , S ( n ) = ( S n , S ) , respectively. Further, it is not hard to compute that every subnetworksuccessively admits a reaction vector balanced equilibrium (1 , (cid:62) , which is consistent withProperty . Finally, according to Remark and Theorem , the equilibrium x ∗ is locallyasymptotically stable.
6. Conclusion.
From the above discussion, the conclusions can be reached that it is possi-ble to capture the stability properties of some large-scale CRNs through the proposed networkdecomposition technique. Some sufficient conditions are presented to suggest stability if thenetwork can be decomposed into a complex balanced subnetwork and a few 1-dimensionalsubnetworks (or a few two-species subnetworks) whether the species are shared among sub-networks or not. The results are finally applied to autocatalytic networks for validation.
Appendix.Appendix 1: CRNs and related stability resutlts.
In this section, some basic conceptsassociated with CRNs are reviewed. For those special CRNs with stability property, theavailable Lyapunov functions are recalled.
A.1. CRNs.
Consider a network consisting of n species S , · · · , S n and r chemical reac-tions. We denote the i th ( i = 1 , · · · , r ) reaction as(A.1) n (cid:88) j =1 v ji S j → n (cid:88) j =1 v (cid:48) ji S j , where v ji , v (cid:48) ji ∈ Z ≥ are the complexes of reactant and product, respectively.We give some definitions related to CRNs [12] that will be used throughout the paper. Definition A.1. (CRN) . A CRN is composed of three finite sets: a set of species S = { S , · · · , S n } ; a set of complexes C = (cid:83) ri =1 { v · i , v (cid:48)· i } and the j th entry of v · i is the stoichiometriccoefficient of S j in this complex, based on which (A.1) is simply written as v .i → v (cid:48) .i ; a set of reactions R = { v · → v (cid:48)· , · · · , v · r → v (cid:48)· r } , satisfying ∀ v · i ∈ C , v · i → v · i / ∈ R but ∃ v (cid:48)· i , s.t. v · i → v (cid:48)· i ∈ R or v (cid:48)· i → v · i ∈ R .The triple ( S , C , R ) or N (cid:44) ( S , C , R ) is usually used to indicate a CRN. Definition A.2. (weakly reversible and reversible CRN)
A CRN ( S , C , R ) is weakly re-versible if for each reaction v · i → v (cid:48)· i ∈ R , there exists a chain of reactions, which startsfrom v (cid:48)· i and ends with v · i , i.e., v (cid:48)· i → v · i ∈ R , v · i → v · i ∈ R , · · · , v · i m → v · i ∈ R . Inparticular, if ∀ v · i → v (cid:48)· i ∈ R , there is v (cid:48)· i → v · i ∈ R , we say that this network is reversible. Definition A.3. (stoichiometric subspace) . For a CRN ( S , C , R ) , the linear subspace S (cid:44) span { v (cid:48)· − v · , · · · , v (cid:48)· r − v · r } is called the stoichiometric subspace of the network, and itsdimension dim S is called the dimension of the network. Definition A.4. (stoichiometric compatibility class).
Given a CRN ( S , C , R ) , for x ∈ R n ≥ ,we say the sets S ( x ) (cid:44) { x + ξ | ξ ∈ S } , ¯ S + ( x ) (cid:44) S ( x ) (cid:84) R n ≥ and S + ( x ) (cid:44) This manuscript is for review purposes only. S ( x ) (cid:84) R n> are the stoichiometric compatibility class, nonnegative and positive stoichiomet-ric compatibility class of x , respectively. Definition A.5. (deficiency)
For a CRN ( S , C , R ) , δ = |C| − l − dim S is defined as the deficiency of the network with |C| and l represent the numbers of complexesand of linkage classes in the network, respectively. When a CRN obeys the mass action law, the reaction rate may be measured according tothe power law with respect to the species concentrations, e.g., for the i th reaction v · i → v (cid:48)· i the reaction rate is Ξ i ( x ) = k i x v · i (cid:44) k i d (cid:89) j =1 x v ji j , where k i ∈ R > is the reaction rate constant, and x ∈ R n ≥ represents the concentration vectorof the species S i with each element x i . Definition A.6. (MAS).
A CRN ( S , C , R ) with mass-action kinetics is called an MAS, rep-resented by the quadruple M (cid:44) ( S , C , R , K ) , where K = { k , · · · , k r } is the set of reaction rateconstants. In the context, we are mainly concerned with MASs. Denote the stoichiometric matrix of N = ( S , C , R ) by Γ ∈ Z n × r with the i th column Γ · i = v (cid:48)· i − v · i , termed reaction vector, andthe r -dimensional non-negative vector function of reaction rate by Ξ( x ) with each componentΞ i ( x ). Then the dynamics of M = ( S , C , R , K ) follows(A.2) d x d t = ΓΞ( x ) , x ∈ R n ≥ . Definition A.7. (equilibrium)
For an M = ( S , C , R , K ) , a positive point x ∗ ∈ R n> is anequilibrium in M if it satisfies ΓΞ( x ∗ ) = 0 . An MAS that admits an equilibrium is a balancedMAS. Following the notion of balanced CRN, we introduce a new concept, termed as generalizedbalancing, which first emerged in [15] in a stochastic version. In this paper, we extended thisnotion to the deterministic case.
Definition A.8. (generalized balancing)
We say an M = ( S , C , R , K ) is generalized bal-anced at x ∗ ∈ R n> if there exists a set of tuples of subsets of reaction set R , denoted as { ( L i , R i ) i ∈ A } with (cid:91) i ∈ A L i = (cid:91) i ∈ A R i = R such that ∀ i ∈ A (cid:88) v · i → v (cid:48)· i ∈ L i k i ( x ∗ ) v · i = (cid:88) v · i → v (cid:48)· i ∈ R i k i ( x ∗ ) v · i (A.3) holds true. This manuscript is for review purposes only. RN DECOMPOSITION FOR STABILITY 15
Remark
A.9. Definition A.8 declares that generalized balancing encompasses the followingcases, • detailed balancing with the set A induced by every reaction, i.e., A = { i | ( v · i −→ v (cid:48)· i , v (cid:48)· i −→ v · i ) v · i −→ v (cid:48)· i ∈R } , L i = (cid:8) v · i → v (cid:48)· i ∈ R (cid:9) , R i = (cid:8) v (cid:48)· i → v · i ∈ R (cid:9) , • complex balancing with the set A induced by any z ∈ C , i.e., A = { z | z ∈ C} , L z = (cid:8) v · i → v (cid:48)· i ∈ R | v · i = z (cid:9) ,R z = (cid:8) v · i → v (cid:48)· i ∈ R | v (cid:48)· i = z (cid:9) , • reaction vector balancing [5] with the set A induced by every given reaction vector η ∈ Z n , i.e., A = { η | η := v (cid:48)· i − v · i , ∀ i } , L η = (cid:8) v · i → v (cid:48)· i ∈ R | v (cid:48)· i − v · i = η (cid:9) ,R η = (cid:8) v · i → v (cid:48)· i ∈ R | v (cid:48)· i − v · i = − η (cid:9) . • any combinatorial forms of the above three types of balancing and other possibilities.Note that a generalized balanced MAS must be a balanced MAS. The former two classes ofbalancing are often encountered in the literature with detailed balancing request the structureof the network to be reversible while complex balanced network to be weakly reversible. It isobvious that a detailed balanced network must be a complex balanced one, and not vice versa.Let us consider the following example to illustrate the concept of reaction vector balancing. Example A.1.
Consider the network describing the autophosphorylation of Aurora B ki-nase, given by E k (cid:47) EP k (cid:111) , E + EP k (cid:47) (cid:47) EP . (A.4) where both trans- and cis-autophosphorylation reactions happens, and
E, EP are proteins.There are complexes, linkage classes, and dim S = 1 , which implies that the deficiency δ = 4 − − . Let S = E, S = EP , there are two reaction vectors ( − , (cid:62) , (1 , − (cid:62) inthe network. If this system has a positive reaction vector balanced equilibrium ( x ∗ , x ∗ ) , thenit should fulfill the following equality k x ∗ + k x ∗ x ∗ = k x ∗ , i.e., ( x ∗ , x ∗ ) = ( c, ck k − ck ) with c < k k a positive constant. A.2. Some known Lyapunov functions.
As stated in the Introduction, for some spe-cial CRNs, such as those with weakly reversible structure and those with dimension 1, theirstability (under some moderate conditions) has been well studied and the corresponding Lya-punov functions are also proposed. For the former under complex balancing condition, thewell-known pseudo-Helmholtz free energy function [18] is an available one, which takes G ( x ) = n (cid:88) j =1 (cid:18) x ∗ j − x j − x j ln x ∗ j x j (cid:19) , x ∈ R n> . (A.5) This manuscript is for review purposes only.
For any 1-dimensional M = ( S , C , R , K ), an available Lyapunov function [8] is given by f ( x ) = (cid:90) γ ( x )0 ln ˜ u ( y † ( x ) + αω )d α, x ∈ R n> (A.6)with ˜ u satisfy h ( x, ˜ u ) = 0 and h ( x, u ) defined by h ( x, u ) = (cid:88) { i | β i > } ( k i x v · i ) (cid:18) β i − (cid:88) j =0 u j (cid:19) + (cid:88) { i | β i < } ( k i x v · i ) (cid:18) − − (cid:88) j = β i u j (cid:19) , (A.7)where ω ∈ R n \{ n } is a set of bases of S , and β i ∈ Z \{ } satisfies v (cid:48)· i − v · i = β i ω, i = 1 , · · · , r ,and moreover, γ ∈ C ( R n> ; R > ), y † ∈ C ( R n> ; R n> ) are constrained by x = y † ( x ) + γ ( x ) ω and γ ( x + bω ) = γ ( x ) + b , ∀ b ∈ R , respectively. Appendix 2: Detailed proofs.
In the appendix, we provide detailed proofs for some resultsappearing in section 3, section 4 and section 5.
B. The proof of Theorem 3.
Denote the state of M ( p ) by x ( p ) ∈ R n p > . Since ∀ p (cid:54) = q ∈ { , , ..., (cid:96) } , E p := S (0) (cid:84) S ( p ) (cid:54) = ∅ , from Definition 1, we can get the state of M : x = (cid:18) (cid:78) S i ∈S (0) x (0) i , (cid:78) (cid:96)p =1 (cid:78) S j ∈S ( p ) \ E p x ( p ) j (cid:19) (cid:62) where x i , x j represent the species S i , S j . Mean-while, according to the decomposition, the dynamics of the M can be represented as˙ x = (cid:96) (cid:88) p =0 r p (cid:88) i =1 k ( p ) i x v ( p ) · i (cid:18) v (cid:48) ( p ) · i − v ( p ) · i (cid:19) , (B.1)where r p is the number of reactions in M ( p ) .Since every M ( p ) for p = 1 , · · · , (cid:96) is reaction vector balanced, plus condition (2) listed inTheorem 3, without loss of generality, the reactions in M ( p ) are supposed to be classified intotwo groups, which respectively are denoted by L ω p = { l : v (cid:48) ( p ) · l − v ( p ) · l = ( (cid:78) S i ∈ E p , ω (cid:62) p ) (cid:62) } and R ω p = { l : v (cid:48) ( p ) · l − v ( p ) · l = − ( (cid:78) S i ∈ E p , ω (cid:62) p ) (cid:62) } . Since M ( p ) is 1-dimensional, we can construct aLyapunov candidate for M ( p ) which is a special form of (A.6), f ( p ) ( x ( p ) ) = (cid:88) S i ∈ E p (cid:18) x ∗ i ( p ) − x ( p ) i − x i ln x ∗ i ( p ) x ( p ) i (cid:19) + (cid:90) γ (˜ x ( p ) )0 ln ˜ u p ( y † (˜ x ( p ) ) + ω p t ) dt, where ˜ x ( p ) = ( (cid:78) S j ∈S ( p ) \ E p x ( p ) j ) (cid:62) , represents the components that are not involved in thecomplex balanced M (0) , γ ∈ C ( R n p −| E p | > ; R > ) , y † ∈ C ( R n p −| E p | > ; R n p −| E p | > ) and ω p sharethe same meanings with (A.6).Combing the Lyapunov function (A.5) for complex balanced M (0) , which is representedas f (0) ( x (0) ) = (cid:88) S i ∈S (0) (cid:18) x ∗ i (0) − x (0) i − x (0) i ln x ∗ i (0) x (0) i (cid:19) , (B.2) This manuscript is for review purposes only.
RN DECOMPOSITION FOR STABILITY 17 it is easy to derive the following function f ( x ) = f (0) ( x (0) ) + (cid:96) (cid:88) p =1 (cid:90) γ (˜ x ( p ) )0 ln ˜ u p ( y † (˜ x ( p ) ) + ω p t ) dt. (B.3)The next step is to verify that f ( x ) is able to be a Lyapunov function for M . Becauseof the continuity of ˜ u p (˜ x ( p ) ) and (3), we know that, for p = 1 , · · · , (cid:96) , there must exist aneighborhood of ˜ x ∗ ( p ) , denoted by N (˜ x ∗ ( p ) ), s.t. ∀ ˜ x ( p ) ∈ N (˜ x ∗ ( p ) ), it holds ω (cid:62) p ∇ ˜ u p (˜ x ( p ) ) > . Since ˜ u p (˜ x ( p ) ) = exp { ω (cid:62) p ∇ f ( p ) (˜ x ( p ) ) } , we have ∇ ˜ u p (˜ x ( p ) ) = ˜ u p ∇ f ( p ) (˜ x ( p ) ) ω p . From Property 1, x ∗ is a generalized equilibrium point in M , then ∀ x ∈ D ( x ∗ ) with D ( x ∗ ) (cid:44) S + ( x ∗ ) (cid:92) (cid:0) R n > (cid:96) (cid:79) p =1 N (˜ x ∗ ( p ) ) (cid:1) , and ∀ µ ∈ S , there is µ (cid:62) ∇ f ( x ) µ = µ (0) (cid:62) diag( (cid:79) S i ∈S (0) /x ∗ (0) i ) µ (0) (cid:124) (cid:123)(cid:122) (cid:125) g ≥ + (cid:96) (cid:88) p =1 µ ( p ) (cid:62) ∇ f ( p ) (˜ x ( p ) ) µ ( p ) = g + (cid:96) (cid:88) p =1 µ ( p ) (cid:62) µ ( p ) ω (cid:62) p ω p ω (cid:62) p ∇ f ( p ) (˜ x ( p ) ) ω p = g + (cid:96) (cid:88) p =1 µ ( p ) (cid:62) µ ( p ) ω (cid:62) p ω p ω (cid:62) p ∇ ˜ u p (˜ x ( p ) )˜ u p ≥ , with the equality holding if and only if µ = n . Thus f ( x ) is strictly convex in this region.Because of the fact ∇ f ( x ∗ ) = (cid:18)(cid:0) Ln x (0) x ∗ (0) (cid:1) (cid:62) , (cid:78) (cid:96)p =1 ω (cid:62) p ω (cid:62) p ω p ln ˜ u (˜ x ( p ) ) (cid:19) (cid:62) | x = x ∗ = (cid:62) n , we obtain f ( x ) ≥ f ( x ∗ ) = 0 for x ∈ D ( x ∗ ).Then we need to prove ˙ f ( x ) ≤
0. Since ∂f ( x ) ∂x ( p ) = ∂f ( p ) ( x ( p ) ) ∂x ( p ) , p = 1 , · · · , (cid:96) , we have˙ f ( x ) = ∇ (cid:62) f ( x ) ˙ x = (cid:18) Ln x (0) x ∗ (0) (cid:19) (cid:62) · r (cid:88) i =1 k (0) i x v (0) · i (cid:18) v (cid:48) (0) · i − v (0) · i (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) ˙ f (0) ( x (0) ) + (cid:96) (cid:88) p =1 r p (cid:88) i =1 k ( p ) i x v ( p ) · i β i ω (cid:62) p ∇ f ( p ) (cid:124) (cid:123)(cid:122) (cid:125) ˙ f ( p ) ( x ( p ) ) , This manuscript is for review purposes only. where the second part satisfies˙ f ( p ) ( x ( p ) ) ≤ r p (cid:88) i =1 k ( p ) i x v ( p ) · i (exp { β i ω (cid:62) p ∇ f ( p ) ( x ( p ) ) } −
1) = 0with the equality holding only if ω (cid:62) p ∇ f ( p ) ( x ( p ) ) = 0, that is (cid:81) S i ∈ E p x ( p ) i x ∗ i ( p ) ˜ u (˜ x ( p ) ) = 1, whichmeans x ( p ) = x ∗ ( p ) . The remaining proof is similar to that of Theorem 2 and we omit it here. (cid:3) C. The proof of Theorem 8.
As we know, the complex balanced M (0) has a Lyapunovfunction in the form of f (0) ( x (0) ) = (cid:88) S i ∈S (0) ( x ∗ i − x i − x i ln x ∗ i x i )to render the asymptotic stability of x (0) . Besides, ∀ p = 1 , · · · , (cid:96) , there is a candidate Lya-punov function for M ( p ) , which is given by f ( p ) ( x ( p ) ) = ( x ∗ i − x i − x i ln x ∗ i x i ) + (cid:90) x j x ∗ j w − jp ln t b ( p ) c ( p ) ij (cid:80) L ωp k ( p ) l t v ( p ) jl dt, where S i ∈ E p , S j ∈ S ( p ) \ E p , and c ( p ) ij = x ∗ i a ( p ) (cid:80) Rωp k ( p ) l x ∗ v ( p ) ili = x ∗ j b ( p ) (cid:80) Lωp k ( p ) l x ∗ v ( p ) jlj .Integrating the above information, we can construct the following function f ( x ) = f (0) ( x (0) ) + (cid:96) (cid:88) p =1 (cid:90) x j x ∗ j w − jp ln t b ( p ) c ( p ) ij (cid:80) L ωp k ( p ) l t v ( p ) jl dt. (C.1)Ulteriorly, we aim to state f ( x ) is localy positive definite. By computing the Hessianmatrix of f ( x ), we derive ∇ f ( x ) = diag (cid:18) (cid:79) S i ∈S (0) x ∗ − i , (cid:79) S j ∈S\S (0) w − jp (cid:80) L ωp k ( p ) l ( b ( p ) − v ( p ) jl ) x v jl ( p ) − j (cid:80) L ωp k ( p ) l x v ( p ) jl j (cid:19) , where (cid:78) stands for the Cartesian product. Since w − jp (cid:80) L ωp k ( p ) l ( b ( p ) − v ( p ) jl ) x v ( p ) jl − j is contin-uous, by the third condition in Theorem 8, there is a neighborhood of x ∗ j for S j ∈ S\S (0) ,denoted by N ( x ∗ j ), such that ∀ x j ∈ N ( x ∗ j ), we have w − jp (cid:80) L ωp k ( p ) l ( b ( p ) − v ( p ) jl ) x ∗ v ( p ) jl − j > . Thus, ∀ x ∈ (cid:0) R n > (cid:78) S j ∈S\S (0) N ( x ∗ j ) (cid:1) (cid:84) S + ( x ∗ ), we get ∇ f ( x ) >
0. Moreover, it is easy toobtain ∇ f ( x ∗ ) = (cid:78) S i ∈S (0) ln x i x ∗ i , (cid:78) S j ∈S\S (0) w − jp ln x ∗ j b ( p ) c ( p ) ij (cid:80) Lωp k ( p ) l x ∗ j v ( p ) jl (cid:62) = (cid:62) n , which im-plies f ( x ) is lower bounded by f ( x ∗ ) = 0. Then we continue to certify ˙ f ( x ) ≤
0. Here, we
This manuscript is for review purposes only.
RN DECOMPOSITION FOR STABILITY 19 use the similar form of (B.1) to represent the dynamics of M although the parameters andcomplexes different, and derive˙ f ( x ) = (cid:18) (cid:79) S i ∈S (0) ln x i x ∗ i (cid:19) (cid:62) r (cid:88) i =1 k (0) i x v (0) · i (cid:18) v (cid:48) (0) · i − v (0) · i (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) ˙ f (0) ( x (0) ) + (cid:96) (cid:88) p =1 r p (cid:88) i =1 k ( p ) i x v ( p ) · i (cid:18) v (cid:48) ( p ) · i − v ( p ) · i (cid:19) (cid:62) ∇ f ( p ) (cid:124) (cid:123)(cid:122) (cid:125) ˙ f ( p ) ( x ( p ) ) . It is obvious that ˙ f (0) ( x (0) ) ≤ x i = x ∗ i for S i ∈ S (0) ,since M (0) is complex balanced. Further, since every M ( p ) , p = 1 , · · · , (cid:96) , is reaction vectorbalanced, according to Lemma 4, there is ˙ f ( p ) ( x ( p ) ) ≤ x ( p ) = x ∗ ( p ) , which directly follows (9). Consequently it states that ˙ f ( x ) ≤ f ( x ) = 0 ifand only if x = x ∗ . (cid:3) D. The proof of Property 2.
In terms of the reactions pair ( R i,j , R j,i ), the dynamics ofthe considered MAS can be represented by˙ x i = (cid:88) S j ∈S (cid:18) (cid:88) R j,i k j,α i x j x α i − i − (cid:88) R i,j k i,α j x i x α j − j (cid:19) (D.1)where i = 1 , · · · , n . Then we explain the sufficiency. Assume every M ( p ) has a reaction vectorbalanced equilibirum x ∗ ( p ) = ( x ∗ i , x ∗ j ) ∈ R > , for p = 1 , · · · , (cid:96) , i, j ∈ { , · · · , n } . Based oncondition (3) in Definition 5, we might as well assume i is fixed and j = 1 , · · · , m , thus, thereis (1) c j k ,α i = k j,α i , ∀ j = 1 , · · · , m i , for all reaction rates in R ji , (2) c j k i,α = k i,α j , ∀ j =1 , · · · , m i , for all reaction rates in R ij . So, Eq. (D.1) can be transformed into˙ x i = m i (cid:88) j =1 c j (cid:18) (cid:88) R j,i k ,α i x j x α i − i − (cid:88) R i,j k i,α x i x α j − j (cid:19) , (D.2)which implies each ( R i,j , R j,i ) has a reaction vector balanced equilibrium ( x ∗ i , x ∗ j ) , j = 1 , · · · , m .Namely, we can know that for ( R i,j , R j,i ) with different j s and fixed i , they share the sameconcentration x ∗ i at their respective equilibria. As a result, the compounded concentration x ∗ = ( x ∗ , · · · , x ∗ n ) from all x ∗ ( p ) is a reaction vector balanced equilibirum of M . In reverse,the necessity is obvious. (cid:3) REFERENCES [1]
M. A. Alradhawi and D. Angeli , New approach to the stability of chemical reaction networks: Piece-wise linear in rates lyapunov functions , IEEE T. Automat. Contr., 61 (2016), pp. 76–89.[2]
D. F. Anderson, G. Craciun, M. Gopalkrishnan, and C. Wiuf , Lyapunov functions, stationary dis-tributions, and non-equilibrium potential for reaction networks , Bull. Math. Biol., 77 (2015), pp. 1744–1767.[3]
D. Angeli , A tutorial on chemical reaction network dynamics , Eur. J. Control., 15 (2009), pp. 398–406.[4]
U. Barenholz, D. Davidi, E. Reznik, Y. M. Baron, N. Antonovsky, E. Noor, and R. Milo , Design principles of autocatalytic cycles constrain enzyme kinetics and force low substrate saturationat flux branch points , eLife., 6 (2017).
This manuscript is for review purposes only. [5]
D. Cappelletti and B. Joshi , Graphically balanced equilibria and stationary measures of reactionnetworks. , SIAM J. Appl. Dyn. Syst., 17 (2018), pp. 2246–2175.[6]
G. Craciun and M. Feinberg , Multiple equilibria in complex chemical reaction networks: Ii. the species-reaction graph , SIAM J. Appl. Math., 66 (2006), pp. 1321–1338.[7]
J. Decraene, G. G. Mitchell, and B. Mcmullin , Crosstalk and the cooperation of collectively auto-catalytic reaction networks , IEEE Congr. Evol. Comput.(CEC), (2009), pp. 2249–2256.[8]
Z. Fang and C. Gao , Lyapunov function partial differential equations for chemical reaction networks:Some special cases , SIAM J. Appl. Dyn. Syst., 18 (2019), pp. 1163–1199.[9]
M. Feinberg , Complex balancing in general kinetic systems , Arch. Ration. Mech. Anal., 49 (1972),pp. 187–194.[10]
M. Feinberg , Chemical reaction network structure and the stability of complex isothermal reactors-i. thedeficiency zero and deficiency one theorems , Chem. Eng. Sci., 42 (1987), pp. 2229–2268.[11]
M. Feinberg , Chemical reaction network structure and the stability of complex isothermal reactors-ii.multiple steady states for networks of deficiency one , Chem. Eng. Sci., 43 (1988), pp. 1–25.[12]
M. Feinberg , The existence and uniqueness of steady states for a class of chemical reaction networks ,Arch. Ration. Mech. Anal., 132 (1995), pp. 311–370.[13]
E. Gross, H. A. Harrington, N. Meshkat, and A. Shiu , Joining and decomposing reaction networks. ,J. Math. Biol., (2020), pp. 1–49.[14]
N. Gruning, H. Lehrach, and M. Ralser , Regulatory crosstalk of the metabolic network , Trends.Biochem. Sci., 35 (2010), pp. 220–227.[15]
L. Hoessly , Stationary distributions via decomposition of stochastic reaction networks , arXiv: Probabil-ity, (2019).[16]
W. Hordijk, J. Hein, and M. Steel , Autocatalytic sets and the origin of life , Entropy, 12 (2010),pp. 1733–1742.[17]
W. Hordijk and M. Steel , Detecting autocatalytic, self-sustaining sets in chemical reaction systems. ,J. Theoret. Biol., 227 (2004), pp. 451–461.[18]
F. Horn and R. Jackson , General mass action kinetics , Arch. Ration. Mech. Anal., 47 (1972), pp. 81–116.[19]
M. D. Johnston, S. David, and S. Gabor , Dynamical equivalence and linear conjugacy of chemicalreaction networks: New results and methods , MATCH Commun. Math. Comput. Chem., 68 (2012),pp. 443–468.[20]
M. D. Johnston and D. Siegel , Linear conjugacy of chemical reaction networks , J. Math. Chem., 49(2011), pp. 1263–1282.[21]
S. A. Kauffman , At home in the universe: The search for laws of self-organization and complexity ,Leonardo., 29 (1995).[22]
M. Ke, Z. Fang, and C. Gao , Complex balancing reconstructed to the asymptotic stability of mass-actionchemical reaction networks with conservation laws , SIAM J. Appl. Math., 79 (2019), pp. 55–74.[23]
Y. Khaluf, C. Pinciroli, G. Valentini, and H. Hamann , The impact of agent density on scalability incollective systems: noise-induced versus majority-based bistability , Swarm. Intell., 11 (2017), pp. 155–179.[24]
C. Lopezotin and T. Hunter , The regulatory crosstalk between kinases and proteases in cancer , Nat.Rev. Cancer, 10 (2010), pp. 278–292.[25]
Y. Lu and C. Gao , Lyapunov function pdes method to the stability of some chemical reaction networks. ,arXiv:2005.01146, (2020).[26]
R. Plasson, A. Brandenburg, L. Jullien, and H. Bersini , Autocatalyses , J. Phys. Chem. A., 115(2011), p. 8073.[27]
G. Szederk´enyi and K. M. Hangos , Finding complex balanced and detailed balanced realizations ofchemical reaction networks , J. Math. Chem., 49 (2011), pp. 1163–1179.[28]
A. van der Schaft, S. Rao, and B. Jayawardhana , On the mathematical structure of balanced chemi-cal reaction networks governed by mass action kinetics , SIAM J. Appl. Math., 73 (2013), pp. 953–973.[29]
N. Virgo, T. Ikegami, and S. Mcgregor , Complex autocatalysis in simple chemistries , Artif. Life.,22 (2016), pp. 1–15.[30]
S. Wu, Y. Lu, and C. Gao , Lyapunov function partial differential equations for stability analysis of aclass of chemical reaction networks. , in 21st IFAC World Congress in Berlin, to appear, 2020., in 21st IFAC World Congress in Berlin, to appear, 2020.