A GraphTheoretical Approach for the Analysis and Model Reduction of ComplexBalanced Chemical Reaction Networks
Abstract
In this paper we derive a compact mathematical formulation describing the dynamics of chemical reaction networks that are complexbalanced 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 complexbalanced networks and for deriving stability properties of such networks. We propose a method for model reduction of complexbalanced 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, zerodeficiency networks, persistence conjecture, equilibria, Schur complement.
1 Introduction
The analysis and control of largescale chemical reaction networks poses a main challenge to bioengineering and systems biology. Chemical reaction networks involving several hundreds of chemical species and reactions are common in living cells [22]. 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 lefthand (substrate) and righthand (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 massaction kinetics chemical reaction networks, which is characterized by the assumption of the existence of an equilibrium for the reaction rates; a socalled 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 [29] is the fact that the dynamics of a detailedbalanced 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 [24] to bear close similarity with consensus algorithms for symmetric multiagent systems. (In fact, it is shown in [29] that the socalled complexaffinities asymptotically reach consensus.) Furthermore, as shown in [17], the framework can be readily extended from mass action kinetics to (reversible) MichaelisMenten 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 wellknown 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 [29] by considering the substantially larger class of complexbalanced reaction networks. A chemical reaction network is called complexbalanced 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 complexbalanced networks was first introduced in [16] and studied in detail in [11, 15, 7, 25, 9]. These systems have also been called as toric dynamical systems in the literature (see [7]).
An example of a complexbalanced network is the model of Tcell interactions due to [19] (see also [27]) depicted in Figure 1.
This chemical reaction network model arises in immunology and was proposed by McKeithan in order to explain the selectivity of Tcell interactions. With reference to Figure 1, and represent a Tcell receptor and a peptidemajor 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 Tcell 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:
(1) 
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 complexbalanced. On the other hand, all reactions in this network are irreversible, and thus McKeithan’s network is not detailedbalanced.
The main aim of this paper is to show how the compact mathematical formulation of detailedbalanced chemical reaction networks derived in [29] can be extended to complexbalanced networks (such as McKeithan’s network). Indeed, the crucial difference between detailedbalanced and complexbalanced 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 multiagent dynamics). In particular, complexbalanced 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 complexbalanced networks share important common characteristics with those of detailedbalanced 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 complexstoichiometric matrix is injective as in the McKeithan’s network case (see Remark 4.7), the same asymptotic stability results were obtained before in [27]. Furthermore, while for detailedbalanced networks it has been shown in [29] that all equilibria are in fact thermodynamic equilibria in this paper the similar result will be proved that all equilibria of a complexbalanced network are complexequilibria. Similar results have already been proved in [16]; however, the proofs presented in the current paper are much more concise and insightful as compared to those presented in [16].
Furthermore, based on our formulation of complexbalanced networks exhibiting a balanced weighted Laplacian matrix associated to the graph of complexes, we will propose a technique for modelreduction of complexbalanced networks. This technique is similar to the Kron reduction method for model reduction of resistive electrical networks described in [18]; 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 largescale 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 complexequilibria and complexbalanced networks and then derive our formulation. In Section 4, we derive equilibrium and stability properties of complexbalanced networks using our formulation. In Section 5, we propose a model reduction method for complexbalanced 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 blockdiagonal 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 timederivative 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 elementwise product and the vector as the elementwise 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 complexbalanced 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]).
2.1 Stoichiometry
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 rowvector 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 [16], 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 complexstoichiometric matrix of the network. Note that by definition all elements of the matrix are nonnegative 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. [5], 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
(3) 
with the incidence matrix of the graph of complexes.
3 The dynamics of complexbalanced 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 complexbalanced 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
(4) 
where is the th element of the complexstoichiometric 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 complexstoichiometric 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
(5) 
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
(6) 
A similar expression of the dynamics corresponding to mass action kinetics, in less explicit form, was already obtained in [27].
3.2 Complexbalanced networks
We now define a class of reaction networks known as complexbalanced networks. This class was first defined in the work of Horn and Jackson (see p. 92 of [16]). We first define a complexequilibrium of a reaction network.
Definition 3.1.
Consider a chemical reaction network with dynamics given by the equation 3. A vector of concentrations is called a complexequilibrium if . Furthermore, a chemical reaction network is called complexbalanced if it admits a complexequilibrium.
It is easy to see that any complexequilibrium 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 complexequilibrium. Observe that
(7) 
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 complexequilibrium, 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 complexequilibrium, every complex involved in the network is at equilibrium.
In [15] and [9], conditions have been derived for a chemical reaction network governed by mass action kinetics to be complexbalanced.
Remark 3.2.
A thermodynamically balanced or detailedbalanced 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. [29]. Such networks are necessarily reversible. Clearly every thermodynamically balanced network is complexbalanced.
We now rewrite the dynamical equations for complexbalanced networks governed by mass action kinetics in terms of a known complexequilibrium. It will be shown that such a form of equations has advantages in deriving stability properties of and also a model reduction method for complexbalanced networks. Recall equation (6) for general mass action reaction networks:
(8) 
Assume that the network is complexbalanced with a complexequilibrium . Define
where denotes the column of . Equation (8) can be rewritten as
(9) 
where . Note that since , also . Furthermore since is a complexequilibrium, we have
Hereafter, we refer to as the weighted Laplacian of the graph of complexes associated with the given complexbalanced 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 complexstoichiometric 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
(10) 
Remark 3.3.
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 complexbalanced networks
In this section, we make use of the compact mathematical formulation (9) in order to derive properties of equilibria and stability of complexbalanced networks.
4.1 Equilibria
Our first result is a characterization of the set of all equilibria of a complexbalanced network in terms of a known equilibrium.
Theorem 4.1.
Consider a complexbalanced network governed by mass action kinetics. Let denote the stoichiometric matrix and assume that is a complexequilibrium for the network. The following hold:

is another equilibrium for the network iff .

Every equilibrium of the network is a complexequilibrium.
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.
Lemma 4.2.
Let be a balanced weighted Laplacian matrix as before. Then for any
Moreover
Proof.
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
(11) 
for all , with equality if and only if . Hence
(12)  
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 .
Proof.
(of Theorem 4.1) The dynamics of the complexbalanced 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 .
(1.)
(Only If): Assume that is an equilibrium, that is
(13) 
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 viceversa. This in turn implies that for every linkage class , consists of equal entries, or in other words it can be written as
where .
Since is a complexequilibrium, . This implies that for , . Now, by evaluating the RHS of (10) at , we have
Remark 4.3.
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 zerodeficiency networks presented in [13]. In the next subsection, we define zerodeficiency networks and prove that every zerodeficiency network that admits an equilibrium is complexbalanced. 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 [13] 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 complexequilibrium . This dependency turns out to be very minor, strengthening the importance of this matrix for the analysis of the network. Indeed, consider any other complexequilibrium . Then , or equivalently . Hence for the th connected component of the complex graph we have , or equivalently, since ,
(14) 
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 Zerodeficiency networks
We now introduce the notion of zerodeficient chemical reaction networks mentioned in Remark 4.3. This notion was introduced in the work of Feinberg [11] in order to relate the stoichiometry of a given network to the structure of the associated graph of complexes.
Definition 4.4.
The deficiency of a chemical reaction network with complexstoichiometric matrix , incidence matrix and stoichiometric matrix is defined as
(15) 
A reaction network has zerodeficiency if .
Note that zerodeficiency is equivalent with ker im, or with the mapping
being injective.
Remark 4.5.
The deficiency of a chemical reaction network has been defined in a different way in [11]. Denote by the number of linkage classes of a given chemical reaction network. Note that as explained in Section 3.3. In [11], deficiency is defined as
(16) 
It is easy to see that definitions (15) and (16) are equivalent.
We now prove that every zerodeficiency network that admits an equilibrium is complexbalanced. Consequently all the results that we state for a complexbalanced network also hold for a zerodeficiency network that admits an equilibrium.
Lemma 4.6.
If a chemical reaction network is zerodeficient and admits an equilibrium, then it is complexbalanced.
Proof.
Consider a zerodeficient network with complexstoichiometric matrix and incidence matrix . Let denote an equilibrium for the given network. Then and hence by zerodeficiency . Consequently is a complexequilibrium and the network is complexbalanced.
The above lemma has been stated and proved earlier in [11, Theorem 4.1, p. 192] in a different and more lengthy manner.
Remark 4.7.
For McKeithan’s network it is easily seen that itself is already injective, thus implying zerodeficiency.
4.3 Asymptotic stability
We now show global asymptotic stability of complexbalanced networks.
Theorem 4.8.
Consider a complexbalanced 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
(17) 
Proof.
Define
(18) 
Observe that
(19) 
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,
(20) 
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
(21) 
and
(22) 
Observe 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 nonzero, and therefore strictly negative, offdiagonal 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).
Remark 4.9.
Remark 4.10.
By proving that is forward invariant with respect to (9), we have actually proved the persistence conjecture stated in [1, p. 1488] of complexbalanced networks obeying the law of conservation of mass. Persistence conjecture roughly corresponds to the requirement of nonextinction of species concentration of a complexbalanced network. Forward invariance of with respect to (6) has been proved in [27, Section VII] for the case when the complexstoichiometric matrix is injective. Proving invariance of is an intricate problem for general mass action kinetics; see e.g. [4] 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 complexbalanced chemical reaction network, there exists a unique complexequilibrium in defined by equation (17). The proof that we provide for this result is very similar to the proof of the zerodeficiency theorem provided in [13] and is based on the following proposition therein. Recall from the Introduction that the product is defined elementwise.
Proposition 4.11.
Let be a linear subspace of , and let . Then there is a unique element , such that .
Proof.
See proof of [13, Proposition B.1, pp. 361363].
Theorem 4.12.
Proof.
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, modelorder reduction is still underdeveloped. The singular perturbation method and quasi steadystate 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 [14], the QSSA approach is extended by considering higherorder approximation in the computation of quasi steadystate. Sunnåker et al. in [26] 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 [23], a method to compute the upperbound of the error is proposed based on sumofsquares programming. The application of these techniques to general kinetics laws, such as MichaelisMenten, poses a significant computational problem.
In this section, we propose a novel and simple method for model reduction of complexbalanced 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 [18]; see also [28]. Moreover, the resulting reducedorder 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.
Proposition 5.1.
Consider a complexbalanced network with a complexequilibrium 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 offdiagonal elements are nonnegative.

and , where .