A Graph-Theoretical Approach for the Analysis and Model Reduction of Complex-Balanced Chemical Reaction Networks
In this paper we derive a compact mathematical formulation describing the dynamics of chemical reaction networks that are complex-balanced and are governed by mass action kinetics. The formulation is based on the graph of (substrate and product) complexes and the stoichiometric information of these complexes, and crucially uses a balanced weighted Laplacian matrix.
It is shown that this formulation leads to elegant methods for characterizing the space of all equilibria for complex-balanced networks and for 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 the computation 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-engineering and systems biology. Chemical reaction networks involving several hundreds of chemical species and reactions are common in living cells . The complexity of their dynamics is further increased by the fact that the chemical reaction rates are intrinsically nonlinear. In particular mass action kinetics, which is the most basic form of expressing the reaction rates, corresponds to differential equations which are polynomial in the concentrations of the chemical species. Chemical reaction network dynamics thus are a prime example of complex networked dynamical systems, and there is a clear need to develop a mathematical framework for handling their complexity.
One of the issues in formalizing complex chemical reaction network dynamics is the fact that their graph representation is not immediate; due to the fact that chemical reactions usually involve more than one substrate chemical species and more than one product chemical species. This problem is resolved by associating the complexes of the reactions, i.e. the left-hand (substrate) and right-hand (product) sides of each reaction, with the vertices of a graph, and the reactions with the edges
Using the graph of complexes formalism, we have developed in [29, 24] a compact mathematical formulation for a class of mass-action kinetics chemical reaction networks, which is characterized by the assumption of the existence of an equilibrium for the reaction rates; a so-called thermodynamic equilibrium. This corresponds to the thermodynamically justified assumption of microscopic reversibility, with the resulting conditions 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 feature of the formulation of  is the fact that the dynamics of a detailed-balanced chemical reaction network is completely specified by a symmetric weighted Laplacian matrix, defined by the graph of complexes and the equilibrium constants, together with an energy function, which is subsequently used for the stability analysis of the network. In particular, the resulting dynamics is shown  to bear close similarity with consensus algorithms for symmetric multi-agent systems. (In fact, it is shown in  that the so-called complex-affinities asymptotically reach consensus.) Furthermore, as shown in , the framework can be readily 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 reaction network models, including the McKeithan network to be explained shortly afterwards. Motivated by such examples we will extend in this paper the results of  by considering the substantially larger class of complex-balanced reaction networks. A chemical reaction network is called complex-balanced if there exists a vector of species concentrations at which the combined rate of outgoing reactions from any complex is equal to the combined rate of incoming reactions to the complex, or in other words each of the complexes involved in the network is at equilibrium. The notion of complex-balanced networks was first introduced in  and studied in detail in [11, 15, 7, 25, 9]. These systems have also been called as toric dynamical systems in the literature (see ).
This chemical reaction network model arises in immunology and was proposed by McKeithan in order to explain the selectivity of T-cell interactions. With reference to Figure 1, and represent a T-cell receptor and a peptide-major histocompatibility complex (MHC) respectively and is a complex for the network. For , represent various intermediate complexes in the phosphorylation and other intermediate modifications of the T-cell receptor ; represents the rate constant of the step of the phosphorylation and is the dissociation rate of the complex. In the following, we denote by the concentration of a species participating in a chemical reaction network. The governing law of the reaction network is the law of mass action kinetics. This leads to the following set of differential equations describing the rate of change of concentrations of various species involved in the network:
Observe that if the left hand side of each of the above equations is set to zero, all the concentrations for can be parametrized in terms of . Since all the rate constants and dissociation constants are positive, it is easy to see that there exists a set of positive concentrations for which 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 network is not detailed-balanced.
The main aim of this paper is to show how the compact mathematical formulation of detailed-balanced chemical reaction networks derived in  can be extended to complex-balanced networks (such as McKeithan’s network). Indeed, the crucial difference between detailed-balanced and complex-balanced networks will turn out to be that in the latter case the Laplacian matrix is not symmetric anymore, but still balanced (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 consensus dynamics with balanced Laplacian matrix. Exploiting this formulation it will be shown how the dynamics of complex-balanced networks share important common characteristics with those of detailed-balanced networks, including a similar characterization of the set of all equilibria and the same stability result stating that the system converges to an equilibrium point uniquely determined by the initial condition. For the particular case when the complex-stoichiometric matrix is injective as in the McKeithan’s network case (see Remark 4.7), the same asymptotic stability results were obtained before in . Furthermore, while for detailed-balanced networks it has been shown in  that all equilibria are in fact thermodynamic equilibria in this paper the similar result will be proved that all equilibria of a complex-balanced network are complex-equilibria. Similar results have already been proved in ; however, the proofs presented in the current paper are much more concise and insightful as compared to those presented in .
Furthermore, based on our formulation of complex-balanced networks exhibiting a balanced weighted Laplacian matrix associated to the graph of complexes, we will propose a technique for model-reduction of complex-balanced networks. This technique is similar to the Kron reduction method for model reduction of resistive electrical networks described in ; see also [10, 28]. Our technique works by deleting complexes from the graph of complexes associated with the network. In other words, our reduced network has fewer complexes and usually fewer reactions as compared to the original network, and yet the behavior of a number of significant metabolites in the reduced network is approximately the same as in the original network. Thus our model reduction method is useful from a computational point of view, specially when we need to deal with models of large-scale chemical reaction networks. Mathematically our approach is based on the result that the Schur complement (with respect to the deleted complexes) of the balanced weighted Laplacian matrix of the full graph of complexes is again a balanced weighted Laplacian matrix corresponding to the reduced graph of complexes.
The paper is organized as follows. In Section 2, we introduce tools from stoichiometry of reactions and graph theory that are required to derive our mathematical formulation. In Section 3, we explain mass action kinetics, define complex-equilibria and complex-balanced networks and then derive our formulation. In Section 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 presents conclusions based on our results.
Notation: The space of dimensional real vectors is denoted by , and the space of real matrices by . The space of dimensional real vectors consisting of all strictly positive entries is denoted by and the space of dimensional real vectors consisting of all nonnegative entries is denoted by . The rank of a real matrix is denoted by . dim denotes the dimension of a set . Given , denotes the diagonal matrix with diagonal entries ; this notation is extended to the block-diagonal case when are real square matrices. Furthermore, and denote the kernel and span respectively of a real matrix . If denotes a linear subspace of , then denotes its orthogonal subspace (with respect to the standard Euclidian inner product). denotes a vector of dimension with all entries equal to 1. The time-derivative of a vector depending on time will be usually denoted by .
Define the mapping as the mapping whose -th component is given as Similarly, define the mapping as the mapping whose -th component is given as Also, define for any vectors the vector as the element-wise product and the vector as the element-wise quotient . Note that with these notations and .
2 Chemical reaction network structure
In this section, we introduce the tools necessary in order to derive our mathematical formulation of the dynamics of complex-balanced networks. First we introduce the concept of stoichiometric matrix of a reaction network. We then define the concept of a graph of complexes, which was first introduced in the work of Horn & Jackson and Feinberg ([16, 15, 12, 13]).
Consider a chemical reaction network involving chemical species (metabolites), among which chemical reactions take place. The basic structure underlying the dynamics of the vector of concentrations of the chemical species is given by the balance laws , where is an matrix, called the stoichiometric matrix. The elements of the vector are commonly called the (reaction) fluxes or rates. The stoichiometric matrix , which consists of (positive and negative) integer elements, captures the basic conservation laws of the reactions. It already contains useful information about the network dynamics, independent of the precise form of the reaction rate . Note that the reaction rate depends 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 an example shown in Figure 2.
Note that the above reaction has 3 irreversible and 1 reversible reaction leading to in total 5 reactions among four species , , and . Since the stoichiometric matrix maps the space of reactions to the space of species, it has dimension . The entry of corresponding to the row and column is obtained by subtracting the number of moles of the species on the product side from that on the substrate side for the reaction. Thus
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 -dimensional row-vector such that , then the quantity is a conserved quantity or a conserved moiety for the dynamics for all possible reaction rates . Indeed, . Later on, in Remark 3.3, it will be shown that law of conservation of mass leads to a conserved moiety of a chemical reaction network.
2.2 Graph of complexes
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 , introducing the space of complexes. The set of complexes of a chemical reaction network is simply defined as the union of all the different left- and righthand sides (substrates and products) of the reactions in the network. Thus, the complexes corresponding to the network (2) are , , and .
The expression of the complexes in terms of the chemical species is formalized by an matrix , whose -th column captures the expression of the -th complex in the chemical species. For example, for the network depicted in Figure 2,
We will call the complex-stoichiometric matrix of the network. Note that by definition all elements of the matrix are non-negative integers.
Since the complexes are the left- and righthand sides of the reactions, they can be naturally associated with the vertices of a directed graph with edges corresponding to the reactions. Formally, the reaction between the -th (reactant) and the -th (product) complex defines a directed edge with tail vertex being the -th complex and head vertex being the -th complex. The resulting graph will be called the graph of complexes.
Recall, see e.g. , that any graph is defined by its incidence matrix . This is a matrix, being the number of vertices and being the number of edges, with -th element equal to if vertex is the tail vertex of edge and if vertex is the head vertex of edge , while otherwise.
Obviously, there is a close relation between the matrix and the stoichiometric matrix . In fact, it is easily checked that
with the incidence matrix of the graph of complexes.
3 The dynamics of complex-balanced networks governed by mass action kinetics
In this section, we first recall the dynamics of species concentrations of reactions governed by mass action kinetics. We then define complex-balanced networks and derive a compact mathematical formulation for their dynamics.
3.1 The general form of mass action kinetics
Recall that for a chemical reaction network, the relation between the reaction rates and species concentrations depends on the governing laws of the reactions involved in the network. In this section, we explain this relation for reaction networks governed by mass action kinetics. In other words, if denotes the vector of reaction rates and denotes the species concentration vector, we show how to construct . The reaction rate of the -th reaction of a mass action chemical reaction network, from a substrate complex to a product complex , is given as
where is the -th element of the complex-stoichiometric matrix , and is the rate constant of the -th reaction. Without loss of generality we will assume throughout that for every , the constant is positive.
This can be rewritten in the following way. Let denote the column of the complex-stoichiometric matrix corresponding to the substrate complex of the -th reaction. Using the mapping as defined at the end of the Introduction, the mass action reaction equation (4) for the -th reaction from substrate complex to product complex can be rewritten as
Based on the formulation in (5), we can describe the complete reaction network dynamics as follows. Let the mass action rate for the complete set of reactions be given by the vector . For every , define
and . Thus if there is no reaction , then . Define the weighted adjacency matrix of the graph of complexes as the matrix with -th element , where . Furthermore, define , where is the diagonal matrix whose -th element is equal to the sum of the elements of the -th column of . Let denote the incidence matrix of the graph of complexes associated with the network. By definition of , we have . It can be verified that the vector for the mass action reaction rate vector is equal to , where the mapping has been defined at the end of the Introduction. Hence the dynamics can be compactly written as
A similar expression of the dynamics corresponding to mass action kinetics, in less explicit form, was already obtained in .
3.2 Complex-balanced networks
We now define a class of reaction networks known as complex-balanced networks. This class was first defined in the work of Horn and Jackson (see p. 92 of ). We first define a complex-equilibrium of a reaction network.
Consider a chemical reaction network with dynamics given by the equation 3. A vector of concentrations is called a complex-equilibrium if . Furthermore, a chemical reaction network 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 way round need not be true (since need not be injective). We now explain the physical interpretation of a complex-equilibrium. Observe that
consists of equations where denotes the number of complexes. Among all the reactions that the complex is involved in, let denote the set of all the reactions for which is the substrate complex and let denote the set of all reactions for which is the product complex. The of equations (7) can now be written as
It follows that at a complex-equilibrium, the combined rate of outgoing reactions from any complex is equal to the combined rate of incoming reactions to the complex. In other words, at a complex-equilibrium, every complex involved in the network is at equilibrium.
A thermodynamically balanced or detailed-balanced chemical reaction network is one for which there exists a vector of positive species concentrations at which each of the reactions of the network is at equilibrium, that is, , see e.g. . Such networks are necessarily reversible. Clearly every thermodynamically balanced network is complex-balanced.
We now rewrite the dynamical equations for complex-balanced networks governed by mass action kinetics in terms of a known complex-equilibrium. It will be shown that such a form of equations has advantages in deriving stability properties of and also a model reduction method for complex-balanced networks. Recall equation (6) for general mass action reaction networks:
Assume that the network is complex-balanced with a complex-equilibrium . Define
where denotes the column of . Equation (8) can be rewritten as
where . Note that since , also . Furthermore since is a complex-equilibrium, we have
Hereafter, we refer to as the weighted Laplacian of the graph of complexes associated with the given complex-balanced network. Both the row and column sums of the weighted Laplacian are equal to zero
3.3 The linkage classes of a graph of complexes
A linkage class of a chemical reaction network is a maximal set of complexes such that is connected by reactions to for every . It can be easily verified that the number of linkage classes of a network, which is equal to the number of connected components of the graph of complexes corresponding to the network, is given by rank (the number of linkage classes in the terminology of [16, 12, 13]). The graph of complexes is connected, i.e., there is one linkage class in the network if and only if ker.
Assume that the reaction network has linkage classes. Assume that the linkage class has reactions between complexes. Partition , and matrices according to the various linkage classes present in the network as follows:
where for , , and denote the complex-stoichiometric matrix, incidence matrix and the weighted Laplacian matrices corresponding to equilibrium concentration respectively for the linkage class. Let denote the stoichiometric matrix of the linkage class. It is easy to see that . Observe that equation (9) can be written as
The law of conservation of mass states that there exists , such that for . This implies that is a conserved moiety for the dynamics , for all forms of the reaction rate . Indeed, implies , since .
4 Equilibria and stability of complex-balanced networks
In this section, we make use of the compact mathematical formulation (9) in order to derive properties of equilibria 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 a known equilibrium.
Consider a complex-balanced network governed by mass action kinetics. Let denote the stoichiometric matrix and assume that is a complex-equilibrium for the network. The following hold:
is another equilibrium for the network iff .
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 for the proof of Theorem 4.8.
Let be a balanced weighted Laplacian matrix as before. Then for any
Let denote the element of and let denote the negative of the element of corresponding to the row and column. Note that . Using the expression can be rewritten as
Furthermore, since the exponential function is strictly convex
for all , with equality if and only if . Hence
since is balanced.
Furthermore, equality occurs in inequality (12) only when each of the terms within the summation on the left hand side is equal to the corresponding term within the summation on the right hand side. Since if the complex reacts to the complex, it follows that for each such , which is equivalent to .
(of Theorem 4.1) The dynamics of the complex-balanced network with complexes are given by
where and are as defined in the previous section. Let denote the incidence matrix of the graph of complexes associated with the network. Assume that the reaction network has linkage classes. Assume that the linkage class has reactions between complexes. Partition , and matrices according to the various linkage classes present in the network as in Section 3.3. Define for .
(Only If): Assume that is an equilibrium, that is
Define . Premultiplying equation (13) with , we get
Hence by Lemma 4.2, and thus .
(If) Assume that . Hence for every linkage class , , or, equivalently, . This implies that if the complex reacts to the complex or vice-versa. This in turn implies that for every linkage class , consists of equal entries, or in other words it can be written as
Since is a complex-equilibrium, . This implies that for , . Now, by evaluating the RHS of (10) at , we have
The steps followed in the proof of the above theorem are very similar to the proof of the characterization of the space of equilibria of a class of networks known as zero-deficiency networks presented in . In the next subsection, we define zero-deficiency networks and prove that every zero-deficiency network that admits an equilibrium is complex-balanced. We emphasize here that the proof of Theorem 4.1 that is presented in this paper is much more simple as compared to similar proofs provided in  due to the use of the properties of the weighted Laplacian in the present manuscript.
One may wonder to what extent the balanced weighted Laplacian matrix depends on the choice of the complex-equilibrium . This dependency turns out to be very minor, strengthening the importance of this matrix for the analysis of the network. Indeed, consider any other complex-equilibrium . Then , or equivalently . Hence for the -th connected component of the complex graph we have , or equivalently, since ,
for some constant . Thus from the definition of , it follows that for some positive constant . Hence, on every connected component of the graph of complexes, the balanced weighted Laplacian matrix is unique up to multiplication by a positive constant.
4.2 Zero-deficiency networks
We now introduce the notion of zero-deficient chemical reaction networks mentioned in Remark 4.3. This notion was introduced in the work of Feinberg  in order to relate the stoichiometry of a given network to the structure of the associated graph of complexes.
The deficiency of a chemical reaction network with complex-stoichiometric matrix , incidence matrix and stoichiometric matrix is defined as
A reaction network has zero-deficiency if .
Note that zero-deficiency is equivalent with ker im, or with the mapping
The deficiency of a chemical reaction network has been defined in a different way in . Denote by the number of linkage classes of a given chemical reaction network. Note that as explained in Section 3.3. In , deficiency is defined as
We now prove that every zero-deficiency network that admits an equilibrium is complex-balanced. Consequently all the results that we state for a complex-balanced network also hold for a zero-deficiency network that admits an equilibrium.
If a chemical reaction network is zero-deficient and admits an equilibrium, then it is complex-balanced.
Consider a zero-deficient network with complex-stoichiometric matrix and incidence matrix . Let denote an equilibrium for the given network. Then and hence by zero-deficiency . Consequently is a complex-equilibrium and the network is complex-balanced.
The above lemma has been stated and proved earlier in [11, Theorem 4.1, p. 192] in a different and more lengthy manner.
For McKeithan’s network it is easily seen that itself is already injective, thus implying zero-deficiency.
4.3 Asymptotic stability
We now show global asymptotic stability of complex-balanced networks.
Consider a complex-balanced network with stoichiometric matrix , an equilibrium and dynamics given by equation 9. Assume that the network obeys the law of conservation of mass stated in Remark 3.3. Then for every initial concentration vector , the species concentration vector converges as to an element of the set
and is proper, i.e., for every real the set is compact. Let and denote the elements of and respectively. From the strict concavity of the logarithmic function,
with equality occuring only when . Putting in equation (20), we get
This implies that
with equality occuring only when , thus proving (19). Properness of is readily checked.
We first prove that
Defining we thus obtain
and the statement follows from Lemma 4.2.
We now prove that is forward invariant with respect to (9). Assume by contradiction that this is not the case, and that for some and . Without loss of generality assume that the species with concentration is present in at least one complex which is involved in a reaction. Let be the subset of complexes which contain , and denote by the matrix containing the columns of corresponding to the complexes in . Hence all elements of the -th row are different from zero. Then it follows that if and thus
where is the -th row vector of . This inequality is due to the fact that the terms corresponding to the positive -th diagonal element of the weighted Laplacian matrix are all zero, while there is at least one term corresponding to a non-zero, and therefore strictly negative, off-diagonal element of . It is easy to see that the inequality also holds at points arbitrarily close to the boundaries of the positive orthant except the origin. Recall that according to the law of conservation of mass, there exists such that is a conserved moiety. Consequently, starting from an initial concentration vector , the state trajectory can not reach the origin. This implies that must be forward invariant with respect to (9).
By proving that 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 conservation of mass. Persistence conjecture roughly corresponds to the requirement of nonextinction of species concentration of a complex-balanced network. Forward invariance of with respect to (6) has been proved in [27, Section VII] for the case when the complex-stoichiometric matrix is injective. Proving invariance of is an intricate problem for general mass action kinetics; see e.g.  and the references quoted therein.
4.4 Equilibrium concentration corresponding to an initial concentration
In this section, we show that corresponding to every positive stoichiometric compatibility class (see equation (2) for a definition) of a complex-balanced chemical reaction network, there exists a unique complex-equilibrium in defined by equation (17). The proof that we provide for this result is very similar to the proof of the zero-deficiency theorem provided in  and is based on the following proposition therein. Recall from the Introduction that the product is defined element-wise.
Let be a linear subspace of , and let . Then there is a unique element , such that .
See proof of [13, Proposition B.1, pp. 361-363].
With reference to Proposition 4.11, define , and observe that . By Proposition 4.11, there exists a unique such that . Define . It follows that , i.e., . Also, since , which is an invariant set of the dynamics. Together with Theorem 4.8 it follows that the state trajectory with initial condition converges to the equilibrium .
5 Model reduction
For biochemical reaction networks, model-order reduction is still underdeveloped. The singular perturbation method 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ärdin , the QSSA approach is extended by considering higher-order approximation in the computation of quasi steady-state. Sunnåker et al. in  proposed a reduction method by identifying variables that can be lumped together and can be used to infer back the original state. In Prescott & Papachristodoulou , a method to compute the upper-bound of the error is proposed based on sum-of-squares programming. The application of these techniques to general kinetics laws, such as Michaelis-Menten, poses a significant computational problem.
In this section, we propose a novel and simple method for model reduction of complex-balanced chemical reaction networks governed by mass action kinetics. Our method is based on the Kron reduction method for model reduction of resistive electrical networks described in ; see also . Moreover, the resulting reduced-order model retains the structure of the kinetics and gives result to a reduced biochemical reaction network, which enables a direct biochemical interpretation.
5.1 Description of the method
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.
Consider a complex-balanced network with a complex-equilibrium and weighted Laplacian matrix corresponding to the equilibrium . Let denote the set of vertices of the graph of complexes associated with the network. Then for any subset of vertices , the Schur complement of with respect to the indices corresponding to satisfies the following properties:
All diagonal elements of are positive and off-diagonal elements are nonnegative.
and , where .