A Graph-Theoretical Approach for the Analysis and Model Reduction of Complex-Balanced Chemical Reaction Networks
AA Graph-Theoretical Approach forthe Analysis and Model Reduction ofComplex-Balanced Chemical Reaction Networks
Shodhan Rao ∗ Arjan van der Schaft † Bayu Jayawardhana ‡ December 12, 2012
Abstract
In this paper we derive a compact mathematical formulation describing the dynamics of chemicalreaction networks that are complex-balanced and are governed by mass action kinetics. The formulationis based on the graph of (substrate and product) complexes and the stoichiometric information of thesecomplexes, and crucially uses a balanced weighted Laplacian matrix. It is shown that this formulationleads to elegant methods for characterizing the space of all equilibria for complex-balanced networks andfor deriving stability properties of such networks. We propose a method for model reduction of complex-balanced networks, which is similar to the Kron reduction method for electrical networks and involves thecomputation of Schur complements of the balanced weighted Laplacian matrix.
Keywords:
Weighted Laplacian matrix, linkage classes, zero-deficiency networks, persistence conjecture,equilibria, Schur complement.
The analysis and control of large-scale chemical reaction networks poses a main challenge to bio-engineeringand systems biology. Chemical reaction networks involving several hundreds of chemical species and reactionsare common in living cells [22]. The complexity of their dynamics is further increased by the fact that thechemical reaction rates are intrinsically nonlinear . In particular mass action kinetics, which is the mostbasic form of expressing the reaction rates, corresponds to differential equations which are polynomial inthe concentrations of the chemical species. Chemical reaction network dynamics thus are a prime exampleof complex networked dynamical systems, and there is a clear need to develop a mathematical frameworkfor handling their complexity. ∗ Center for Systems Biology, University of Groningen, email: [email protected] † Johann Bernoulli Institute for Mathematics and Computer Science, University of Groningen, e-mail:
[email protected] ‡ Discrete Technology and Production Automation, University of Groningen, email: [email protected] a r X i v : . [ m a t h . D S ] N ov ne of the issues in formalizing complex chemical reaction network dynamics is the fact that theirgraph representation is not immediate; due to the fact that chemical reactions usually involve more thanone substrate chemical species and more than one product chemical species. This problem is resolved byassociating the complexes of the reactions, i.e. the left-hand (substrate) and right-hand (product) sides ofeach reaction, with the vertices of a graph, and the reactions with the edges . The resulting directed graph,called the graph of complexes , is characterized by its incidence matrix B . Furthermore the stoichiometricmatrix S of the chemical reaction network, expressing the basic balance laws of the reactions, can befactorized as S = ZB , with the complex-stoichiometric matrix Z encoding the expressions of the complexesin the various chemical species. We note that an alternative graph formulation of chemical reaction networksis the species-reaction graph [8, 2, 3, 4], which is a bipartite graph with one part of the vertices correspondingto the species and the remaining part to the reactions, and the edges expressing the involvement of the speciesin the reactions.Using the graph of complexes formalism, we have developed in [29, 24] a compact mathematical formula-tion for a class of mass-action kinetics chemical reaction networks, which is characterized by the assumptionof the existence of an equilibrium for the reaction rates; a so-called thermodynamic equilibrium. Thiscorresponds to the thermodynamically justified assumption of microscopic reversibility, with the resultingconditions on the parameters of the mass action kinetics usually referred to as the Wegscheider conditions.The resulting class of mass action kinetics reaction networks are called (detailed)-balanced . A main featureof the formulation of [29] is the fact that the dynamics of a detailed-balanced chemical reaction networkis completely specified by a symmetric weighted Laplacian matrix , defined by the graph of complexes andthe equilibrium constants, together with an energy function, which is subsequently used for the stabilityanalysis of the network. In particular, the resulting dynamics is shown [24] to bear close similarity with consensus algorithms for symmetric multi-agent systems. (In fact, it is shown in [29] that the so-called complex-affinities asymptotically reach consensus.) Furthermore, as shown in [17], the framework can bereadily extended from mass action kinetics to (reversible) Michaelis-Menten reaction rates.On the other hand, the assumption of existence of a thermodynamical equilibrium requires reversibility of all the reactions of the network, while there are quite a few well-known irreversible chemical reactionnetwork models, including the McKeithan network to be explained shortly afterwards. Motivated by suchexamples we will extend in this paper the results of [29] by considering the substantially larger class of complex-balanced reaction networks. A chemical reaction network is called complex-balanced if there exists avector of species concentrations at which the combined rate of outgoing reactions from any complex is equalto the combined rate of incoming reactions to the complex, or in other words each of the complexes involvedin the network is at equilibrium. The notion of complex-balanced networks was first introduced in [16] andstudied in detail in [11, 15, 7, 25, 9]. These systems have also been called as toric dynamical systems in theliterature (see [7]).An example of a complex-balanced network is the model of T-cell interactions due to [19] (see also [27]) This approach is originating in the work of Horn and Jackson [16]; see also Othmer [21] and [3] for nice expos´es andadditional insights. epicted in Figure 1. This chemical reaction network model arises in immunology and was proposed by T + M (cid:45) C (cid:45) C (cid:45) . . . (cid:45) C i (cid:45) . . . (cid:45) C N k p, k − , k − , k − ,i k − ,N k p, k p, k p,i k p,i +1 k p,N Figure 1: McKeithan’s networkMcKeithan in order to explain the selectivity of T-cell interactions. With reference to Figure 1, T and M represent a T-cell receptor and a peptide-major histocompatibility complex (MHC) respectively and T + M is a complex for the network. For i = 1 , . . . , N , C i represent various intermediate complexes inthe phosphorylation and other intermediate modifications of the T-cell receptor T ; k p,i represents the rateconstant of the i th step of the phosphorylation and k − ,i is the dissociation rate of the i th complex. In thefollowing, we denote by [ A ] the concentration of a species A participating in a chemical reaction network.The governing law of the reaction network is the law of mass action kinetics. This leads to the followingset of differential equations describing the rate of change of concentrations of various species involved in thenetwork: d [ T ] dt = d [ M ] dt = − k p, [ T ][ M ] + N (cid:88) i =0 k − ,i [ C i ] d [ C ] dt = k p, [ T ][ M ] − ( k − , + k p, )[ C ]... d [ C i ] dt = k p,i [ C i − ] − ( k − ,i + k p,i +1 )[ C i ]... d [ C N ] dt = k p,N [ C N − ] − k − ,N [ C N ] (1)Observe that if the left hand side of each of the above equations is set to zero, all the concentrations [ C i ] for i = 0 , . . . , N − , can be parametrized in terms of [ C N ]. Since all the rate constants and dissociation constantsare positive, it is easy to see that there exists a set of positive concentrations { [ T ] , [ M ] , [ C ] , . . . , [ C N ] } forwhich the right hand sides of the equations (1) vanish. This implies that McKeithan’s network is complex-balanced. On the other hand, all reactions in this network are irreversible , and thus McKeithan’s networkis not detailed-balanced.The main aim of this paper is to show how the compact mathematical formulation of detailed-balancedchemical reaction networks derived in [29] can be extended to complex-balanced networks (such as McKei-than’s network). Indeed, the crucial difference between detailed-balanced and complex-balanced networkswill turn out to be that in the latter case the Laplacian matrix is not symmetric anymore, but still bal-anced (in the sense of the terminology used in graph theory and multi-agent dynamics). In particular,complex-balanced chemical networks will be shown to bear close resemblance with asymmetric consensusdynamics with balanced Laplacian matrix. Exploiting this formulation it will be shown how the dynamicsf complex-balanced networks share important common characteristics with those of detailed-balanced net-works, including a similar characterization of the set of all equilibria and the same stability result statingthat the system converges to an equilibrium point uniquely determined by the initial condition. For theparticular case when the complex-stoichiometric matrix Z is injective as in the McKeithan’s network case(see Remark 4.7), the same asymptotic stability results were obtained before in [27]. Furthermore, whilefor detailed-balanced networks it has been shown in [29] that all equilibria are in fact thermodynamic equi-libria in this paper the similar result will be proved that all equilibria of a complex-balanced network arecomplex-equilibria. Similar results have already been proved in [16]; however, the proofs presented in thecurrent paper are much more concise and insightful as compared to those presented in [16].Furthermore, based on our formulation of complex-balanced networks exhibiting a balanced weightedLaplacian matrix associated to the graph of complexes, we will propose a technique for model-reduction ofcomplex-balanced networks. This technique is similar to the Kron reduction method for model reduction ofresistive electrical networks described in [18]; see also [10, 28]. Our technique works by deleting complexesfrom the graph of complexes associated with the network. In other words, our reduced network has fewercomplexes and usually fewer reactions as compared to the original network, and yet the behavior of a numberof significant metabolites in the reduced network is approximately the same as in the original network. Thusour model reduction method is useful from a computational point of view, specially when we need to dealwith models of large-scale chemical reaction networks. Mathematically our approach is based on the resultthat the Schur complement (with respect to the deleted complexes) of the balanced weighted Laplacianmatrix of the full graph of complexes is again a balanced weighted Laplacian matrix corresponding to thereduced graph of complexes.The paper is organized as follows. In Section 2, we introduce tools from stoichiometry of reactions andgraph theory that are required to derive our mathematical formulation. In Section 3, we explain mass actionkinetics, define complex-equilibria and complex-balanced networks and then derive our formulation. InSection 4, we derive equilibrium and stability properties of complex-balanced networks using our formulation.In Section 5, we propose a model reduction method for complex-balanced networks, while Section 6 presentsconclusions based on our results. Notation : The space of n dimensional real vectors is denoted by R n , and the space of m × n real matricesby R m × n . The space of n dimensional real vectors consisting of all strictly positive entries is denoted by R n + and the space of n dimensional real vectors consisting of all nonnegative entries is denoted by ¯ R n + . The rankof a real matrix A is denoted by rank A . dim( V ) denotes the dimension of a set V . Given a , . . . , a n ∈ R ,diag( a , . . . , a n ) denotes the diagonal matrix with diagonal entries a , . . . , a n ; this notation is extended tothe block-diagonal case when a , . . . , a n are real square matrices. Furthermore, ker A and span A denote thekernel and span respectively of a real matrix A . If U denotes a linear subspace of R m , then U ⊥ denotesits orthogonal subspace (with respect to the standard Euclidian inner product). m denotes a vector ofdimension m with all entries equal to 1. The time-derivative dxdt ( t ) of a vector x depending on time t will beusually denoted by ˙ x .efine the mapping Ln : R m + → R m , x (cid:55)→ Ln( x ) , as the mapping whose i -th component is given as(Ln( x )) i := ln( x i ) . Similarly, define the mapping Exp : R m → R m + , x (cid:55)→ Exp( x ) , as the mapping whose i -th component is given as (Exp( x )) i := exp( x i ) . Also, define for any vectors x, z ∈ R m the vector x · z ∈ R m as the element-wise product ( x · z ) i := x i z i , i = 1 , , . . . , m, and the vector xz ∈ R m as the element-wisequotient (cid:0) xz (cid:1) i := x i z i , i = 1 , · · · , m . Note that with these notations Exp( x + z ) = Exp( x ) · Exp( z ) andLn( x · z ) = Ln( x ) + Ln( z ) , Ln (cid:0) xz (cid:1) = Ln( x ) − Ln( z ). In this section, we introduce the tools necessary in order to derive our mathematical formulation of thedynamics of complex-balanced networks. First we introduce the concept of stoichiometric matrix of areaction network. We then define the concept of a graph of complexes, which was first introduced in thework of Horn & Jackson and Feinberg ([16, 15, 12, 13]).
Consider a chemical reaction network involving m chemical species (metabolites), among which r chemicalreactions take place. The basic structure underlying the dynamics of the vector x ∈ R m + of concentrations x i , i = 1 , · · · , m, of the chemical species is given by the balance laws ˙ x = Sv , where S is an m × r matrix,called the stoichiometric matrix . The elements of the vector v ∈ R r are commonly called the (reaction) fluxes or rates . The stoichiometric matrix S , which consists of (positive and negative) integer elements,captures the basic conservation laws of the reactions. It already contains useful information about thenetwork dynamics, independent of the precise form of the reaction rate v ( x ). Note that the reaction ratedepends on the governing law prescribing the dynamics of the reaction network.We now show how to construct the stoichiometric matrix for a reaction network with the help of anexample shown in Figure 2. X + 2 X X (cid:10) (cid:45) X + X (cid:25) X (cid:107) Figure 2: Example of a reaction networkNote that the above reaction has 3 irreversible and 1 reversible reaction leading to in total 5 reactionsamong four species X , X , X and X . Since the stoichiometric matrix maps the space of reactions tothe space of species, it has dimension 4 ×
5. The entry of S corresponding to the i th row and j th columnis obtained by subtracting the number of moles of the i th species on the product side from that on theubstrate side for the j th reaction. Thus S = − − − − − − −
10 0 0 1 1 for the reaction network depicted in Figure 2.In this paper, we focus only on closed chemical reaction networks meaning those without external fluxes.Therefore unless otherwise mentioned, our chemical reaction networks do not have any external fluxes.If there exists an m -dimensional row-vector k such that kS = 0, then the quantity kx is a conservedquantity or a conserved moiety for the dynamics ˙ x = Sv ( x ) for all possible reaction rates v = v ( x ). Indeed, ddt kx = kSv ( x ) = 0. Later on, in Remark 3.3, it will be shown that law of conservation of mass leads to aconserved moiety of a chemical reaction network.Note that for all possible fluxes the solutions of the differential equations ˙ x = Sv ( x ) starting from aninitial state x will always remain within the affine space S x := { x ∈ R m + | x − x ∈ im S } . (2) S x has been referred to as the positive stoichiometric compatibility class (corresponding to x ) in [13, 25, 1]. The structure of a chemical reaction network cannot be directly captured by an ordinary graph. Instead,we will follow an approach originating in the work of Horn and Jackson [16], introducing the space of complexes . The set of complexes of a chemical reaction network is simply defined as the union of all thedifferent left- and righthand sides (substrates and products) of the reactions in the network. Thus, thecomplexes corresponding to the network (2) are X + 2 X , X , 2 X + X and X .The expression of the complexes in terms of the chemical species is formalized by an m × c matrix Z ,whose α -th column captures the expression of the α -th complex in the m chemical species. For example, forthe network depicted in Figure 2, Z = . We will call Z the complex-stoichiometric matrix of the network. Note that by definition all elements of thematrix Z are non-negative integers.Since the complexes are the left- and righthand sides of the reactions, they can be naturally associatedwith the vertices of a directed graph G with edges corresponding to the reactions. Formally, the reaction α −→ β between the α -th (reactant) and the β -th (product) complex defines a directed edge with tail vertexeing the α -th complex and head vertex being the β -th complex. The resulting graph will be called the graph of complexes .Recall, see e.g. [5], that any graph is defined by its incidence matrix B . This is a c × r matrix, c beingthe number of vertices and r being the number of edges, with ( α, j )-th element equal to − α isthe tail vertex of edge j and 1 if vertex α is the head vertex of edge j , while 0 otherwise.Obviously, there is a close relation between the matrix Z and the stoichiometric matrix S . In fact, it iseasily checked that S = ZB, hence ˙ x = ZBv ( x ) (3)with B the incidence matrix of the graph of complexes. In this section, we first recall the dynamics of species concentrations of reactions governed by mass actionkinetics. We then define complex-balanced networks and derive a compact mathematical formulation fortheir dynamics.
Recall that for a chemical reaction network, the relation between the reaction rates and species concentrationsdepends on the governing laws of the reactions involved in the network. In this section, we explain thisrelation for reaction networks governed by mass action kinetics. In other words, if v denotes the vector ofreaction rates and x denotes the species concentration vector, we show how to construct v ( x ). The reactionrate of the j -th reaction of a mass action chemical reaction network, from a substrate complex S j to aproduct complex P j , is given as v j ( x ) = k j m (cid:89) i =1 x Z i S j i , (4)where Z iρ is the ( i, ρ )-th element of the complex-stoichiometric matrix Z , and k j ≥ j -th reaction. Without loss of generality we will assume throughout that for every j , the constant k j is positive.This can be rewritten in the following way. Let Z S j denote the column of the complex-stoichiometricmatrix Z corresponding to the substrate complex of the j -th reaction. Using the mapping Ln : R c + → R c as defined at the end of the Introduction, the mass action reaction equation (4) for the j -th reaction fromsubstrate complex S j to product complex P j can be rewritten as v j ( x ) = k j exp (cid:0) Z T S j Ln( x ) (cid:1) . (5)ased on the formulation in (5), we can describe the complete reaction network dynamics as follows. Letthe mass action rate for the complete set of reactions be given by the vector v ( x ) = (cid:104) v ( x ) · · · v r ( x ) (cid:105) T .For every σ, π ∈ { , . . . , c } , define C πσ := { j ∈ { , . . . , r } | ( σ, π ) = ( S j , P j ) } and a πσ := (cid:80) j ∈C πσ k j . Thus if there is no reaction σ → π , then a πσ = 0. Define the weighted adjacencymatrix A of the graph of complexes as the matrix with ( π, σ )-th element a πσ , where π, σ ∈ { , . . . , c } .Furthermore, define L := ∆ − A , where ∆ is the diagonal matrix whose ( ρ, ρ )-th element is equal to thesum of the elements of the ρ -th column of A . Let B denote the incidence matrix of the graph of complexesassociated with the network. By definition of L , we have Tc L = 0. It can be verified that the vector Bv ( x ) for the mass action reaction rate vector v ( x ) is equal to − L Exp (cid:0) Z T Ln( x ) (cid:1) , where the mappingExp : R c → R c + has been defined at the end of the Introduction. Hence the dynamics can be compactlywritten as˙ x = − ZL Exp (cid:0) Z T Ln( x ) (cid:1) (6)A similar expression of the dynamics corresponding to mass action kinetics, in less explicit form, was alreadyobtained in [27]. We now define a class of reaction networks known as complex-balanced networks. This class was first definedin the work of Horn and Jackson (see p. 92 of [16]). We first define a complex-equilibrium of a reactionnetwork.
Definition 3.1.
Consider a chemical reaction network with dynamics given by the equation ( ) . A vectorof concentrations x ∗ ∈ R m + is called a complex-equilibrium if Bv ( x ∗ ) = 0 . Furthermore, a chemical reactionnetwork is called complex-balanced if it admits a complex-equilibrium. It is easy to see that any complex-equilibrium is an equilibrium for the network, but the other wayround need not be true (since Z need not be injective). We now explain the physical interpretation of acomplex-equilibrium. Observe that Bv ( x ∗ ) = 0 (7)consists of c equations where c denotes the number of complexes. Among all the reactions that the i th complex C i is involved in, let O i denote the set of all the reactions for which C i is the substrate complexand let I i denote the set of all reactions for which C i is the product complex. The i th of equations (7) cannow be written as (cid:88) k ∈I i v k ( x ∗ ) = (cid:88) k ∈O i v k ( x ∗ )t follows that at a complex-equilibrium, the combined rate of outgoing reactions from any complex is equalto the combined rate of incoming reactions to the complex. In other words, at a complex-equilibrium, everycomplex involved in the network is at equilibrium.In [15] and [9], conditions have been derived for a chemical reaction network governed by mass actionkinetics to be complex-balanced. Remark 3.2. A thermodynamically balanced or detailed-balanced chemical reaction network is one for whichthere exists a vector of positive species concentrations x ∗ at which each of the reactions of the network isat equilibrium, that is, v ( x ∗ ) = 0, see e.g. [29]. Such networks are necessarily reversible. Clearly everythermodynamically balanced network is complex-balanced.We now rewrite the dynamical equations for complex-balanced networks governed by mass action kineticsin terms of a known complex-equilibrium. It will be shown that such a form of equations has advantages inderiving stability properties of and also a model reduction method for complex-balanced networks. Recallequation (6) for general mass action reaction networks:˙ x = − ZL Exp (cid:0) Z T Ln( x ) (cid:1) (8)Assume that the network is complex-balanced with a complex-equilibrium x ∗ . Define K ( x ∗ ) := diag ci =1 (cid:0) exp( Z Ti Ln( x ∗ )) (cid:1) where Z i denotes the i th column of Z . Equation (8) can be rewritten as˙ x = − ZLK ( x ∗ ) K ( x ∗ ) − Exp (cid:0) Z T Ln( x ) (cid:1) = − Z L ( x ∗ )Exp (cid:16) Z T Ln (cid:16) xx ∗ (cid:17)(cid:17) (9)where L ( x ∗ ) := LK ( x ∗ ). Note that since Tc L = 0, also Tc L ( x ∗ ) = 0. Furthermore since x ∗ is a complex-equilibrium, we have L ( x ∗ ) c = L ( x ∗ )Exp (cid:16) Z T Ln (cid:16) xx ∗ (cid:17)(cid:17) | x = x ∗ = 0Hereafter, we refer to L ( x ∗ ) as the weighted Laplacian of the graph of complexes associated with the givencomplex-balanced network. Both the row and column sums of the weighted Laplacian L ( x ∗ ) are equal tozero . It is this special property of the weighted Laplacian that we make use of in deriving all the resultsstated further on in this paper. A linkage class of a chemical reaction network is a maximal set of complexes {C , . . . , C k } such that C i isconnected by reactions to C j for every i, j ∈ { , . . . , k } , i (cid:54) = j . It can be easily verified that the numberof linkage classes ( (cid:96) ) of a network, which is equal to the number of connected components of the graph of In the literature on directed graphs (see e.g. [6]), L is called balanced . Note that the matrix L , having zero column sumsbut not zero row sums, is similar to the ‘advection’ set-up considered in [6]. omplexes corresponding to the network, is given by (cid:96) = c − rank( B ) (the number of linkage classes in theterminology of [16, 12, 13]). The graph of complexes is connected, i.e., there is one linkage class in thenetwork if and only if ker( B T ) = span (cid:0) c (cid:1) .Assume that the reaction network has (cid:96) linkage classes. Assume that the i th linkage class has r i reactionsbetween c i complexes. Partition Z , B and L matrices according to the various linkage classes present in thenetwork as follows: Z = (cid:104) Z Z . . . Z (cid:96) (cid:105) B = B . . . B . . . . . . B (cid:96) − . . . . . . B (cid:96) L ( x ∗ ) = diag( L ( x ∗ ) , L ( x ∗ ) , . . . , L (cid:96) − ( x ∗ ) , L (cid:96) ( x ∗ ))where for ( i = 1 , . . . , (cid:96) ), Z i ∈ ¯ R m × c i + , B i ∈ R c i × r i and L i ( x ∗ ) ∈ R c i × c i denote the complex-stoichiometricmatrix, incidence matrix and the weighted Laplacian matrices corresponding to equilibrium concentration x ∗ respectively for the i th linkage class. Let S i denote the stoichiometric matrix of the i th linkage class. Itis easy to see that S i = Z i B i . Observe that equation (9) can be written as˙ x = − (cid:96) (cid:88) i =1 Z i L i ( x ∗ )Exp (cid:16) Z Ti Ln (cid:16) xx ∗ (cid:17)(cid:17) (10) Remark 3.3.
The law of conservation of mass states that there exists u ∈ R m + , such that Z Ti u ∈ span (cid:0) c i (cid:1) for i = 1 , . . . , (cid:96) . This implies that u T x is a conserved moiety for the dynamics ˙ x = ZBv , for all forms of thereaction rate v ( x ). Indeed, Z Ti u ∈ span (cid:0) c i (cid:1) implies u T ZB = 0, since B Ti c i = 0. In this section, we make use of the compact mathematical formulation (9) in order to derive properties ofequilibria and stability of complex-balanced networks.
Our first result is a characterization of the set of all equilibria of a complex-balanced network in terms of aknown equilibrium.
Theorem 4.1.
Consider a complex-balanced network governed by mass action kinetics. Let S ∈ R m × r denote the stoichiometric matrix and assume that x ∗ ∈ R m + is a complex-equilibrium for the network. Thefollowing hold:1. x ∗∗ ∈ R m + is another equilibrium for the network iff S T Ln (cid:0) x ∗∗ x ∗ (cid:1) = 0 .. Every equilibrium of the network is a complex-equilibrium. The proof of this theorem crucially makes use of the following lemma, which will also be the basis forthe proof of Theorem 4.8.
Lemma 4.2.
Let L ( x ∗ ) be a balanced weighted Laplacian matrix as before. Then for any γ ∈ R c γ T L ( x ∗ )Exp( γ ) ≥ Moreover γ T L ( x ∗ )Exp( γ ) = 0 if and only if B T γ = 0 Proof.
Let γ i denote the i th element of γ and let k ij denote the negative of the element of L ( x ∗ ) correspondingto the j th row and i th column. Note that k ij ≥ , i, j = 1 , · · · , c . Using Tc L ( x ∗ ) = 0 the expression − γ T L ( x ∗ )Exp( γ ) can be rewritten as (cid:104)(cid:80) i (cid:54) =1 ( γ i − γ ) k i (cid:80) i (cid:54) =2 ( γ i − γ ) k i . . . (cid:80) i (cid:54) = c ( γ i − γ c ) k ci (cid:105) Exp( γ )= c (cid:88) j =1 (cid:88) i (cid:54) = j ( γ i − γ j ) k ji exp( γ j )Furthermore, since the exponential function is strictly convex( β − α )exp( α ) ≤ exp( α ) − exp( β ) (11)for all α, β , with equality if and only if α = β . Hence − γ T L ( x ∗ )Exp( γ ) = c (cid:88) j =1 (cid:88) i (cid:54) = j ( γ i − γ j ) k ji exp( γ j ) ≤ c (cid:88) j =1 (cid:88) i (cid:54) = j k ji (cid:0) exp( γ i ) − exp( γ j ) (cid:1) (12)= Tc (cid:80) i (cid:54) =1 k i (cid:0) exp( γ i ) − exp( γ ) (cid:1) ... (cid:80) i (cid:54) = c k ci (cid:0) exp( γ i ) − exp( γ c ) (cid:1) = − Tc L ( x ∗ ) T Exp( γ )= − (cid:0) L ( x ∗ ) c (cid:1) T Exp( γ ) = 0 . since L ( x ∗ ) is balanced.Furthermore, equality occurs in inequality (12) only when each of the terms within the summation onthe left hand side is equal to the corresponding term within the summation on the right hand side. Since k ij > , i (cid:54) = j, if the i th complex reacts to the j th complex, it follows that γ i = γ j for each such i, j , which isequivalent to B T γ = 0. (cid:4) Proof. (of Theorem 4.1) The dynamics of the complex-balanced network with c complexes are given by˙ x = − Z L ( x ∗ )Exp (cid:16) Z T Ln (cid:16) xx ∗ (cid:17)(cid:17) here Z ∈ ¯ R m × c + and L ( x ∗ ) ∈ R c × c are as defined in the previous section. Let B ∈ R c × r denote theincidence matrix of the graph of complexes associated with the network. Assume that the reaction networkhas (cid:96) linkage classes. Assume that the i th linkage class has r i reactions between c i complexes. Partition Z , B and L matrices according to the various linkage classes present in the network as in Section 3.3. Define S i := Z i B i for i = 1 , . . . , (cid:96) .( .)( Only If ): Assume that x ∗∗ is an equilibrium, that is − Z L ( x ∗ )Exp (cid:18) Z T Ln (cid:18) x ∗∗ x ∗ (cid:19)(cid:19) = 0 (13)Define γ := Z T Ln (cid:0) x ∗∗ x ∗ (cid:1) . Premultiplying equation (13) with Ln (cid:0) x ∗∗ x ∗ (cid:1) , we get − γ T L ( x ∗ )Exp( γ ) = 0 . Hence by Lemma 4.2, B T γ = 0 and thus S T Ln (cid:0) x ∗∗ x ∗ (cid:1) = 0.( If ) Assume that S T Ln (cid:0) x ∗∗ x ∗ (cid:1) = 0. Hence for every linkage class i = 1 , . . . , (cid:96) , S Ti Ln (cid:0) x ∗∗ x ∗ (cid:1) = 0, or,equivalently, B Ti γ = 0. This implies that γ i = γ j if the i th complex reacts to the j th complex or vice-versa.This in turn implies that for every linkage class i = 1 , . . . , (cid:96) , Z Ti Ln (cid:0) x ∗∗ x ∗ (cid:1) consists of equal entries, or in otherwords it can be written as Z Ti Ln (cid:18) x ∗∗ x ∗ (cid:19) = Γ i c i where Γ i ∈ R .Since x ∗ is a complex-equilibrium, L ( x ∗ ) c = 0. This implies that for i = 1 , . . . , (cid:96) , L i ( x ∗ ) c i = 0. Now,by evaluating the RHS of (10) at x ∗∗ , we have − (cid:96) (cid:88) i =1 Z i L i ( x ∗ )Exp (cid:18) Z Ti Ln (cid:18) x ∗∗ x ∗ (cid:19)(cid:19) = − (cid:96) (cid:88) i =1 exp(Γ i ) Z i L i ( x ∗ ) c i = 0 . ( .) Let x ∗∗ denote an equilibrium as in the proof of the earlier part. Then S T Ln (cid:0) x ∗∗ x ∗ (cid:1) = 0. Weprove that x ∗∗ is a complex-equilibrium. As shown earlier, Z Ti Ln (cid:0) x ∗∗ x ∗ (cid:1) = Γ i c i for i = 1 , . . . , (cid:96) . Since L i ( x ∗ ) c i = 0, this implies that (c.f., the discussion that preceeds (6) and the form in (10)) − Bv ( x ∗∗ ) = L ( x ∗ )Exp (cid:18) Z T Ln (cid:18) x ∗∗ x ∗ (cid:19)(cid:19) = L ( x ∗ )Exp (cid:0) Z T Ln (cid:0) x ∗∗ x ∗ (cid:1)(cid:1) ... L (cid:96) ( x ∗ )Exp (cid:0) Z T(cid:96) Ln (cid:0) x ∗∗ x ∗ (cid:1)(cid:1) = 0From the above equation, it follows that x ∗∗ is a complex-equilibrium. (cid:4) Remark 4.3.
The steps followed in the proof of the above theorem are very similar to the proof of thecharacterization of the space of equilibria of a class of networks known as zero-deficiency networks presentedn [13]. In the next subsection, we define zero-deficiency networks and prove that every zero-deficiencynetwork that admits an equilibrium is complex-balanced. We emphasize here that the proof of Theorem 4.1that is presented in this paper is much more simple as compared to similar proofs provided in [13] due tothe use of the properties of the weighted Laplacian in the present manuscript.One may wonder to what extent the balanced weighted Laplacian matrix L ( x ∗ ) depends on the choiceof the complex-equilibrium x ∗ . This dependency turns out to be very minor, strengthening the importanceof this matrix for the analysis of the network. Indeed, consider any other complex-equilibrium x ∗∗ . Then S T Ln x ∗∗ = S T Ln x ∗ , or equivalently B T Z T Ln x ∗∗ = B T Z T Ln x ∗ . Hence for the i -th connected componentof the complex graph we have B Ti Z Ti Ln x ∗∗ = B Ti Z Ti Ln x ∗ , or equivalently, since ker B Ti = span , Z Ti Ln x ∗∗ = Z Ti Ln x ∗ + c i (14)for some constant c i . Thus from the definition of L , it follows that L i ( x ∗∗ ) = d i L i ( x ∗ ) for some positiveconstant d i . Hence, on every connected component of the graph of complexes, the balanced weightedLaplacian matrix L ( x ∗ ) is unique up to multiplication by a positive constant. We now introduce the notion of zero-deficient chemical reaction networks mentioned in Remark 4.3. Thisnotion was introduced in the work of Feinberg [11] in order to relate the stoichiometry of a given networkto the structure of the associated graph of complexes.
Definition 4.4.
The deficiency δ of a chemical reaction network with complex-stoichiometric matrix Z ,incidence matrix B and stoichiometric matrix S is defined as δ := rank( B ) − rank( ZB ) = rank( B ) − rank( S ) ≥ A reaction network has zero-deficiency if δ = 0 . Note that zero-deficiency is equivalent with ker( Z ) ∩ im( B ) = 0, or with the mapping Z : im B ⊂ R c → R m being injective . Remark 4.5.
The deficiency of a chemical reaction network has been defined in a different way in [11].Denote by (cid:96) the number of linkage classes of a given chemical reaction network. Note that (cid:96) = c − rank( B )as explained in Section 3.3. In [11], deficiency δ is defined as δ := c − (cid:96) − rank( S ) (16)It is easy to see that definitions (15) and (16) are equivalent.e now prove that every zero-deficiency network that admits an equilibrium is complex-balanced. Con-sequently all the results that we state for a complex-balanced network also hold for a zero-deficiency networkthat admits an equilibrium. Lemma 4.6.
If a chemical reaction network is zero-deficient and admits an equilibrium, then it is complex-balanced.Proof.
Consider a zero-deficient network with complex-stoichiometric matrix Z ∈ ¯ R m × c + and incidence matrix B ∈ R c × r . Let x ∗ ∈ R m + denote an equilibrium for the given network. Then Sv ( x ∗ ) = ZBv ( x ∗ ) = 0 andhence by zero-deficiency Bv ( x ∗ ) = 0. Consequently x ∗ is a complex-equilibrium and the network is complex-balanced. (cid:4) The above lemma has been stated and proved earlier in [11, Theorem 4.1, p. 192] in a different and morelengthy manner.
Remark 4.7.
For McKeithan’s network it is easily seen that Z itself is already injective, thus implyingzero-deficiency. We now show global asymptotic stability of complex-balanced networks.
Theorem 4.8.
Consider a complex-balanced network with stoichiometric matrix S ∈ R m × r , an equilibrium x ∗ ∈ R m + and dynamics given by equation ( ) . Assume that the network obeys the law of conservation ofmass stated in Remark 3.3. Then for every initial concentration vector x ∈ R m + , the species concentrationvector x converges as t → ∞ to an element of the set E := { x ∗∗ ∈ R m + | S T Ln( x ∗∗ ) = S T Ln( x ∗ ) } . (17) Proof.
Define G ( x ) = x T Ln (cid:16) xx ∗ (cid:17) + ( x ∗ − x ) T m (18)Observe that G ( x ∗ ) = 0 . We prove that G ( x ) > ∀ x (cid:54) = x ∗ , (19)and is proper , i.e., for every real c > { x ∈ ¯ R m + | G ( x ) ≤ c } is compact. Let x i and x ∗ i denote the i th elements of x and x ∗ respectively. From the strict concavity of the logarithmic function, z − ln( z ) ≥ ∀ z ∈ R + with equality occuring only when z = 1. Putting z = x ∗ i x i in equation (20), we get x ∗ i − x i + x i ln (cid:18) x i x ∗ i (cid:19) ≥ G defined by (18) is a standard Lyapunov function used in chemical reaction network theory (see for example [13]). his implies that G ( x ) = m (cid:88) i =1 (cid:20) x ∗ i − x i + x i ln (cid:18) x i x ∗ i (cid:19)(cid:21) ≥ . with equality occuring only when x = x ∗ , thus proving (19). Properness of G is readily checked.We first prove that˙ G ( x ) := ∂G∂x ( x ) ˙ x ≤ ∀ x ∈ R m + , (21)and ˙ G ( x ) = 0 if and only if x ∈ E . (22)Observe that˙ G ( x ) = − Ln (cid:16) xx ∗ (cid:17) T Z T L ( x ∗ )Exp (cid:16) Z T Ln (cid:16) xx ∗ (cid:17)(cid:17) Defining γ := Z T Ln (cid:0) xx ∗ (cid:1) we thus obtain˙ G ( x ) = − γ T L ( x ∗ )Exp( γ )and the statement follows from Lemma 4.2.We now prove that R m + is forward invariant with respect to (9). Assume by contradiction that this isnot the case, and that x i ( t ) = 0 for some t and i ∈ { , . . . , m } . Without loss of generality assume thatthe species with concentration x i is present in at least one complex which is involved in a reaction. Let C i be the subset of complexes which contain x i , and denote by Z C i the matrix containing the columns of Z corresponding to the complexes in C i . Hence all elements of the i -th row Z C i are different from zero. Thenit follows that (cid:81) mj =1 x Z jC j = 0 if C ∈ C i and thus˙ x i ( t ) = − Z i L ( x ∗ ) (cid:81) mj =1 x j ( t ) Z jC ... (cid:81) mj =1 x j ( t ) Z jCc > , where Z i is the i -th row vector of Z . This inequality is due to the fact that the terms corresponding to thepositive i -th diagonal element of the weighted Laplacian matrix L ( x ∗ ) are all zero, while there is at leastone term corresponding to a non-zero, and therefore strictly negative, off-diagonal element of L ( x ∗ ). It iseasy to see that the inequality also holds at points arbitrarily close to the boundaries of the positive orthant R m + except the origin. Recall that according to the law of conservation of mass, there exists u ∈ R m + suchthat u T x is a conserved moiety. Consequently, starting from an initial concentration vector x ∈ R m + , thestate trajectory x ( · ) can not reach the origin. This implies that R m + must be forward invariant with respectto (9).Since G is proper (in ¯ R m + ) and the state trajectory x ( · ) remains in R m + , (21) implies that x ( · ) is boundedin R m + . Therefore, boundedness of x ( · ), together with equations (21) and (22), imply that the speciesconcentration x converges to an element of the set E by an application of LaSalle’s invariance principle. (cid:4) emark 4.9. The crux of the proofs of Theorems 4.1 and 4.8 is the inequality γ T L ( x ∗ )Exp( γ ) ≥ , for all γ ∈ R . This inequality holds because of balancedness of L and we make use of the convexity of the exponentialfunction in order to prove it. Remark 4.10.
By proving that R m + is forward invariant with respect to (9), we have actually proved the persistence conjecture stated in [1, p. 1488] of complex-balanced networks obeying the law of conservationof mass. Persistence conjecture roughly corresponds to the requirement of nonextinction of species concen-tration of a complex-balanced network. Forward invariance of R m + with respect to (6) has been proved in[27, Section VII] for the case when the complex-stoichiometric matrix Z is injective. Proving invariance of R m + is an intricate problem for general mass action kinetics; see e.g. [4] and the references quoted therein. In this section, we show that corresponding to every positive stoichiometric compatibility class (see equa-tion (2) for a definition) of a complex-balanced chemical reaction network, there exists a unique complex-equilibrium in E defined by equation (17). The proof that we provide for this result is very similar tothe proof of the zero-deficiency theorem provided in [13] and is based on the following proposition therein.Recall from the Introduction that the product x · z ∈ R m is defined element-wise. Proposition 4.11.
Let U be a linear subspace of R m , and let x ∗ , x ∈ R m + . Then there is a unique element µ ∈ U ⊥ , such that (cid:0) x ∗ · Exp( µ ) − x (cid:1) ∈ U .Proof. See proof of [13, Proposition B.1, pp. 361-363]. (cid:4)
Theorem 4.12.
Consider the complex-balanced chemical reaction network with dynamics given by equation ( ) and equilibrium set E . Then for every x ∈ R m + there exists a unique x ∈ E ∩ S x with S x given by ( ) ,and the solution trajectory x ( · ) of ( ) with initial condition x (0) = x converges for t → ∞ to x .Proof. With reference to Proposition 4.11, define U = im S , and observe that U ⊥ = ker S T . By Proposition4.11, there exists a unique µ ∈ ker S T such that x ∗ · Exp( µ ) − x ∈ im S . Define x := x ∗ · Exp( µ ) ∈ R m + . Itfollows that S T µ = S T Ln (cid:0) x x ∗ (cid:1) = 0, i.e., x ∈ E . Also, since x − x ∈ im S , x ∈ S x which is an invariant setof the dynamics. Together with Theorem 4.8 it follows that the state trajectory x ( · ) with initial condition x (0) = x converges to the equilibrium x ∈ E . (cid:4) For biochemical reaction networks, model-order reduction is still underdeveloped. The singular perturbationmethod and quasi steady-state approximation (QSSA) approach are the most commonly used techniques,where the reduced state contains a part of the metabolites of the full model. In the thesis by H¨ardin [14], theQSSA approach is extended by considering higher-order approximation in the computation of quasi steady-state. Sunn˚aker et al. in [26] proposed a reduction method by identifying variables that can be lumpedogether and can be used to infer back the original state. In Prescott & Papachristodoulou [23], a method tocompute the upper-bound of the error is proposed based on sum-of-squares programming. The applicationof these techniques to general kinetics laws, such as Michaelis-Menten, poses a significant computationalproblem.In this section, we propose a novel and simple method for model reduction of complex-balanced chemicalreaction networks governed by mass action kinetics. Our method is based on the Kron reduction methodfor model reduction of resistive electrical networks described in [18]; see also [28]. Moreover, the resultingreduced-order model retains the structure of the kinetics and gives result to a reduced biochemical reactionnetwork, which enables a direct biochemical interpretation.
Our model reduction method is based on reduction of the graph of complexes associated with the network.It is based on the following result regarding Schur complements of weighted Laplacian matrices.
Proposition 5.1.
Consider a complex-balanced network with a complex-equilibrium x ∗ ∈ R m + and weightedLaplacian matrix L ( x ∗ ) ∈ R c × c corresponding to the equilibrium x ∗ . Let V denote the set of vertices ofthe graph of complexes associated with the network. Then for any subset of vertices V r ⊂ V , the Schurcomplement ˆ L ( x ∗ ) of L ( x ∗ ) with respect to the indices corresponding to V r satisfies the following properties:1. All diagonal elements of ˆ L ( x ∗ ) are positive and off-diagonal elements are nonnegative.2. T ˆ c ˆ L ( x ∗ ) = 0 and ˆ L ( x ∗ ) ˆ c = 0 , where ˆ c := c − dim ( V r ) .Proof. ( .) Follows from the proof of [20, Theorem 3.11]; see also [28] for the case of symmetric L .( .) Without loss of generality, assume that the last c − ˆ c rows and columns of L ( x ∗ ) correspond to thevertex set V r . Consider a partition of L ( x ∗ ) given by L ( x ∗ ) = (cid:34) L ( x ∗ ) L ( x ∗ ) L ( x ∗ ) L ( x ∗ ) (cid:35) (23)where L ( x ∗ ) ∈ R ˆ c × ˆ c , L ( x ∗ ) ∈ R ˆ c × ( c − ˆ c ) , L ( x ∗ ) ∈ R ( c − ˆ c ) × ˆ c and L ( x ∗ ) ∈ R ( c − ˆ c ) × ( c − ˆ c ) . By definition,ˆ L ( x ∗ ) = L ( x ∗ ) − L ( x ∗ ) L ( x ∗ ) − L ( x ∗ )Since Tc L ( x ∗ ) = 0, we obtain T ˆ c L ( x ∗ ) + Tc − ˆ c L ( x ∗ ) = 0 T ˆ c L ( x ∗ ) + Tc − ˆ c L ( x ∗ ) = 0Using the above equations, we get T ˆ c ˆ L ( x ∗ ) = T ˆ c (cid:0) L ( x ∗ ) − L ( x ∗ ) L ( x ∗ ) − L ( x ∗ ) (cid:1) = − Tc − ˆ c L ( x ∗ ) + Tc − ˆ c L ( x ∗ ) L ( x ∗ ) − L ( x ∗ ) = 0In a similar way, it can be proved that ˆ L ( x ∗ ) ˆ c = 0. (cid:4) rom the above result, it follows that ˆ L ( x ∗ ) obeys all the properties of the weighted Laplacian of acomplex-balanced chemical network corresponding to the graph of complexes with vertex set V − V r . ThusProposition 5.1 can be directly applied to the graph of complexes, yielding a reduction of the chemicalreaction network by reducing the number of complexes. Consider a complex-balanced reaction networkdescribed in the standard form (9)Σ : ˙ x = − Z L ( x ∗ )Exp (cid:16) Z T Ln (cid:16) xx ∗ (cid:17)(cid:17) Reduction will be performed by deleting certain complexes in the graph of complexes , resulting in a reducedgraph of complexes with weighted Laplacian ˆ L ( x ∗ ). Furthermore, leaving out the corresponding columns ofthe complex-stoichiometric matrix Z one obtains a reduced complex-stoichiometric matrix ˆ Z (with as manycolumns as the remaining number of complexes in the graph of complexes), leading to the reduced reactionnetworkˆΣ : ˙ x = − ˆ Z ˆ L ( x ∗ )Exp (cid:16) ˆ Z T Ln (cid:16) xx ∗ (cid:17)(cid:17) . (24)Note that ˆΣ is again a complex-balanced chemical reaction network governed by mass action kinetics, witha reduced number of complexes and with, in general, a different set of reactions (edges of the graph ofcomplexes). Furthermore, the complex-equilibrium x ∗ of the original reaction network Σ is a complex-equilibrium of the reduced network ˆΣ as well.An interpretation of the reduced network ˆΣ can be given as follows. Consider a subset V r of the set ofall complexes, which we wish to leave out in the reduced network. Consider the partition of L ( x ∗ ) as givenby equation (23) and a corresponding partition of Z given by Z = (cid:104) Z Z (cid:105) , (25)where V r corresponds to the last part of the indices (denoted by 2), in order to write out the dynamics ofΣ as ˙ x = − (cid:104) Z Z (cid:105) (cid:34) L ( x ∗ ) L ( x ∗ ) L ( x ∗ ) L ( x ∗ ) (cid:35) (cid:34) Exp (cid:0) Z T Ln (cid:0) xx ∗ (cid:1)(cid:1) Exp (cid:0) Z T Ln (cid:0) xx ∗ (cid:1)(cid:1)(cid:35) Consider now the auxiliary dynamical system (cid:34) ˙ y ˙ y (cid:35) = − (cid:34) L ( x ∗ ) L ( x ∗ ) L ( x ∗ ) L ( x ∗ ) (cid:35) (cid:34) w w (cid:35) where we impose the constraint ˙ y = 0. It follows that w = −L ( x ∗ ) − L ( x ∗ ) w , leading to the reduced auxiliary dynamics defined by the Schur complement˙ y = − (cid:0) L ( x ∗ ) − L ( x ∗ ) L ( x ∗ ) − L ( x ∗ ) (cid:1) w = − ˆ L ( x ∗ ) w utting back in w = Exp (cid:16) ˆ Z T Ln (cid:0) xx ∗ (cid:1)(cid:17) , making use of ˙ x = Z ˙ y + Z ˙ y = Z ˙ y = ˆ Z ˙ y , we then obtain thereduced network ˆΣ given in (24).We derive the following properties relating Σ and ˆΣ. Proposition 5.2.
Consider the complex-balanced reaction network Σ and its reduced order model ˆΣ givenby ( ) . Denote their sets of equilibria by E , respectively ˆ E . Then E ⊆ ˆ E .Proof. Let B denote the incidence matrix for Σ and let c and ˆ c denote the number of complexes in Σ andˆΣ respectively. Assume that the graph of complexes is connected; otherwise the same argument can berepeated for every component (linkage class). It follows that ker( B T ) = span( c ). Let x ∗∗ ∈ E . We showthat x ∗∗ ∈ ˆ E . Let Z and ˆ Z denote the complex-stoichiometric matrices of Σ and ˆΣ respectively. Let L ( x ∗ )denote the weighted Laplacian of the graph of complexes of Σ corresponding to an equilibrium x ∗ . Let ˆ L ( x ∗ )denote the Schur complement of L ( x ∗ ) corresponding to the reduced model ˆΣ.Since x ∗∗ ∈ E , B T Z T Ln (cid:0) x ∗∗ x ∗ (cid:1) = 0. It follows that Z T Ln (cid:0) x ∗∗ x ∗ (cid:1) ∈ span( c ). This implies that ˆ Z T Ln (cid:0) x ∗∗ x ∗ (cid:1) ∈ span( ˆ c ) since the columns of ˆ Z form a subset of the columns of Z . From Proposition 5.1, it now followsthat ˆ L ( x ∗ )Exp (cid:16) ˆ Z T Ln (cid:0) x ∗∗ x ∗ (cid:1)(cid:17) = 0. Consequently x ∗∗ ∈ ˆ E . This concludes the proof. (cid:4) In this section, we show the effect of our model reduction method on two types of complex-balanced networkswith a single linkage class. In other words, we give an interpretation of our reduced model in terms of itscorresponding full model for two types of networks. Note that deletion of a set of complexes in one linkageclass does not affect the reactions of the other linkage classes of the network.
Type 1:
Full Network: C k (cid:10) k − C k (cid:10) k − C k (cid:10) k − · · · · · · k n − (cid:10) k − ( n − C n (26)Reduced Network: C k k k − k (cid:10) k − k − k − k C k (cid:10) k − · · · · · · k n − (cid:10) k − ( n − C n (27)These are complex-balanced networks with reversible reactions occuring between consecutive elementsof the set of distinct complexes {C , C , . . . , C n } as in (26). The reduced network obtained by deleting thecomplex C can be computed to be (27). The two reactions, C (cid:10) C and C (cid:10) C in the full networkare replaced by one reaction C (cid:10) C in the reduced network. This reaction is again a reversible reactiongoverned by mass action kinetics, with rate constants given by (27).The transient behaviour of the metabolites involved in the complexes of the reduced model will approxi-mately be the same as that of the full model if the metabolites involved in C reach their steady states muchfaster than the remaining metabolites. In this case, we can safely impose the condition that the metabolitesnvolved in C are at constant concentration in order to obtain the reduced model (27) with similar transientbehaviour as that of (26). The rule of induction may be applied in order to further reduce the model bydeleting more complexes.A special case of Type 1 networks is C (cid:10) C . Deletion of the complex C in this case is equivalent todeletion of the linkage class from the network. Such a deletion provides a close approximation to the originalnetwork if the reaction has very little effect on the dynamics of the network, i.e., if the reaction reaches itssteady state much faster than the remaining reactions of the network. C (cid:45) C (cid:45) C (cid:45) C (cid:45) . . . (cid:45) C n k k − k − k − k − n k k k k n Figure 3: Type 2 full network C (cid:45) C (cid:45) C (cid:45) . . . (cid:45) C n k k − k − + k k − k − + k k − nk k k − + k k k n Figure 4: Type 2 reduced network
Type 2:
These are complex-balanced networks between distinct complexes {C , C , . . . , C n } as shown inFigure 3. McKeithan’s network is an example of such networks. The reduced network obtained by deletingcomplex C is shown in Figure 4. Observe that the reaction C → C in the reduced network has a differentrate constant as compared to the full network. The three reactions C → C , C → C and C → C of the fullnetwork are replaced by one reaction C → C in the reduced network. The rate constant for this reaction isgiven in Figure 4. All the remaining reactions of the reduced network occur in the same way as in the fullnetwork.In this case, the transient behaviour of the metabolites involved in the complexes of the reduced modelwill approximately be the same as that of the full model if the metabolites involved in C reach their steadystates much faster than the remaining metabolites. Using the method described in this paper, we can studythe effects of deleting other complexes like C or C n from the model.Observe that for all the types of networks discussed above, it is important to determine which of thecomplexes are to be deleted so that the reduced model approximates the full model reasonably well. Example 5.3.
We have applied the model reduction method described in this paper to the model of T-cellinteractions as in (1). We use the following numerical values: N = 19; k p, = 52; k p, = 49; k p, = 41; k p, = 39; k p, = 37; k p, = 34; k p, = 31; k p, = 29; k p, = 25; k p, = 19; k p, = 16; k p, = 21; k p, = 20; k p, = 19; k p, = 18; k p, = 15; k p, = 24; k p, = 13; k p, = 7; k p, = 5; k − , = 13; k − , = 29; k − , = 0 . k − , = 1 . k − , = 2 . k − , = 2; k − , = 0 . k − , = 0 . k − , = 0 . k − , = 0 . k − , = 0 . k − , = 0 . k − , = 3; k − , = 5; k − , = 1; k − , = 11; k − , = 0 . k − , = 7; k − , = 1; k − , = 17 . [ T ] [ M ] Full ModelReduced Model Full ModelReduced Model
Figure 5: Concentration profiles of T and M The initial value of each of the complexes C i , i = 0 , . . . ,
19 is assumed to be equal to 0.01. The complexes T and M are assumed to have initial concentrations 1 and 2 respectively. We have performed modelreduction by deleting 5 complexes C , C , C , C and C . We have simulated the transient behaviourof the remaining complexes for the first two time units. The simulation results show that there is a goodagreement between the transient behaviour of the concentration of most of such complexes when comparingthe full network to the reduced network. Figure 5 depicts plots of comparison of the concentration profilesof T and M .An interpretation of model reduction of the above example model is as follows. By deleting complexes C , C , C , C and C , we assume that these complexes are at constant concentration. Since deletingthese complexes results in a reduced model that closely mimics the original model, it follows that in the fullmodel, these complexes reach an equilibrium much faster than the remaining complexes, so that assumingthat these complexes are at constant concentrations results in a close approximation of the original model. In this paper, we have provided a compact mathematical formulation for the dynamics of complex-balancednetworks. We have made use of this formulation for the determination of equilibria and the asymptoticstability of such networks. The methods that have been employed are very similar to the ones used in [13],but the difference is that our proofs are much more concise than the ones presented in [13] due to the use ofproperties of balanced weighted Laplacian matrices of complex-balanced networks. By proving the forwardinvariance of the positive orthant with respect to the dynamics of complex-balanced networks, we believethat we have proved the persistence conjecture for complex-balanced networks stated in [1]. Furthermore,we have made use of the formulation in order to derive a model reduction technique for complex-balancednetworks.A main challenge for further research is the extension of our results to chemical reaction networks withexternal fluxes and/or externally controlled concentrations. This will change the stability analysis consider-ably, due to the nonlinearity of the differential equations. Furthermore, it will also lead to scrutinizing theodel reduction technique proposed in this paper from an external (input-output) point of view.
References [1] D.F. Anderson, “A proof of the global attractor conjecture in the single linkage class case”,
SIAM J.Appl. Math. , 71(4), pp. 1487–1508, 2011.[2] D. Angeli, P. De Leenheer, E.D. Sontag, “Graph-theoretic characterizations of monotonicity of chemicalnetworks in reaction coordinates”, J. Math. Biol., 61, pp. 581–616, 2010.[3] D. Angeli, “A tutorial on chemical reaction network dynamics”,
European Journal of Control , 15 (3-4),pp. 398–406, 2009.[4] D. Angeli, P. De Leenheer, E.D. Sontag, “Persistence results for chemical reaction networks with time-dependent kinetics and no global conservation laws”,
SIAM J. Appl. Math. , 71, pp. 128–146, 2011.[5] B. Bollobas,
Modern Graph Theory , Graduate Texts in Mathematics 184, Springer, New York, 1998.[6] A. Chapman and M. Mesbahi, “Advection on graphs”, 50 th IEEE CDC-ECC , Orlando, USA, pp.1461–1466, 2011.[7] G. Craciun, A. Dickenstein, A. Shiu, B. Sturmfels, “Toric dynamical systems”,
Journal of SymbolicComputation , 44, pp. 1551–1565, 2009.[8] G. Craciun, M. Feinberg, “Multiple equilibria in complex chemical reaction networks: I. the injectivityproperty”,
SIAM J. Appl. Math. , 65(5), pp. 1526–1546, 2006.[9] A. Dickenstein and M.P. Mill´an, “How far is complex balancing from detailed balancing”,
Bull. Math.Biol. , 73, pp. 811–828, 2011.[10] F. D¨orfler, F. Bullo, ”Kron Reduction of Graphs with Applications to Electrical Networks”,
IEEETransactions on Circuits and Systems I , 99, pp. 1-14, 2012.[11] M. Feinberg, “Complex balancing in chemical kinetics”,
Arch. Rational Mech. Anal. , 49, pp. 187–194,1972.[12] M. Feinberg, “Chemical reaction network structure and the stability of complex isothermal reactors-I. The deficiency zero and deficiency one theorems”,
Chemical Engineering Science , 43(10), pp. 2229–2268, 1987.[13] M. Feinberg, “The existence and uniqueness of steady states for a class of chemical reaction networks”,
Arch. Rational Mech. Anal. , 132, pp. 311–370, 1995.[14] H.M. H¨ardin,
Handling Biological Complexity: As simple as possible but not simpler , Ph.D. Thesis,Vrije Universiteit Amsterdam, 2010.15] F.J.M. Horn, “Necessary and sufficient conditions for complex balancing in chemical kinetics”,
Arch.Rational Mech. Anal. , 49, pp. 172–186, 1972.[16] F. Horn and R. Jackson, “General mass action kinetics”,
Arch. Rational Mech. Anal. , 47, pp. 81–116,1972.[17] B. Jayawardhana, S. Rao, A. van der Schaft, “Balanced chemical reaction networks governed by generalkinetics”, Proc. 20th Mathematical Theory of Networks and Systems, Melbourne, June 2012.[18] G. Kron,
Tensor Analysis of Networks , Wiley, 1939.[19] T.W. McKeithan, “Kinetic proofreading in T-cell receptor signal transduction”,
Proc. Natl. Acad. Sci.USA , 92, pp. 5042-5046, 1995.[20] N.M.D. Niezink, “Consensus in networked multi-agent systems”, Master’s thesis in Applied Mathe-matics, Faculty of Mathematics and Natural Sciences, University of Groningen, August 2011.[21] H. G. Othmer,
Analysis of Complex Reaction Networks , Lecture Notes, School of Mathematics, Uni-versity of Minnesota, December 9, 2003.[22] B.O. Palsson,
Systems Biology; Properties of Reconstructed Networks , Cambridge University Press,Cambridge 2006.[23] T.P. Prescott, A. Papachristodoulou, “Guaranteed error bounds for structured complexity reductionof biochemical networks”, Journal of Theoretical Biology, vol. 304, pp. 172-182, 2012.[24] S. Rao, B. Jayawardhana, A.J. van der Schaft, “On the graph and systems analysis of reversiblechemical reaction networks with mass action kinetics”, Proc. IEEE American Control Conference,Montreal, June 2012.[25] D. Siegel and D. MacLean, “Global stability of complex balanced mechanisms”,
Journal of Mathemat-ical Chemistry , 27, pp. 89–110, 2000.[26] M. Sunn˚aker, G. Cedersund, M. Jirstrand, “A method for zooming of nonlinear models of biochemicalsystems”, BMC Systems Biology, 5:140, 2011.[27] E.D. Sontag, “Structure and stability of certain chemical networks and applications to the kineticproofreading model of T-cell receptor signal transduction”,
IEEE Trans. Autom. Control , 46(7), pp.1028–1047, 2001.[28] A.J. van der Schaft, “Characterization and partial synthesis of the behavior of resistive circuits at theirterminals”,