Associated with a chemical reaction network is a natural labelled bipartite multigraph termed an SR graph, and its directed version, the DSR graph. These objects are closely related to Petri nets. The construction of SR and DSR graphs for chemical reaction networks is presented. Conclusions about asymptotic behaviour of the associated dynamical systems which can be drawn easily from the graphs are discussed. In particular, theorems on ruling out the possibility of multiple equilibria or stable oscillation in chemical reaction networks based on computations on SR/DSR graphs are presented. These include both published and new results. The power and limitations of such results are illustrated via several examples.
The applicability of such work, particularly in biological contexts, is greatly
increased if only weak assumptions are made about kinetics. Consequently, there
is a growing body of recent work on CRNs with essentially arbitrary kinetics.
It has been shown that examination of Petri nets associated with a CRN allows
conclusions about persistence, that is, whether ω-limit sets of interior points of
Rn≥0 can intersect the boundary of Rn≥0 [
Outline. After some preliminaries, the construction of SR and DSR graphs
is presented, and their relationship to Petri nets is discussed. Some recent results
about multistationarity based on cycle structure in these objects are described.
Subsequently, a new result on monotonicity in CRNs is proved. This result,
Proposition 4, is a graph-theoretic corollary of results in [
Consider the following simple family of CRNs treated in [
A3
B1
B2 2A1
SYS 2 A1 + A2 A2 + A3 A3 + A4
A4
B1 B2
B3 2A1 · · · · · ·
SYS n Ai + Ai+1 Bi,
i = 1, . . . , n + 1 An+2 2A1 The reader may wish to look ahead to Figure 2 to see representations of the SR graphs associated with the first three CRNs in this family. This family will be revisited in Section 7, and the theory to be presented will imply the following conclusions (to be made precise below): when n is even, SYS n does not allow multiple nondegenerate equilibria; when n is odd, SYS n cannot have a nontrivial periodic attractor. Both conclusions require only minimal assumptions about the kinetics. 2.2
In a spatially homogeneous setting, a chemical reaction system in which n reactants participate in m reactions has dynamics governed by the ordinary differential equation x˙ = Γ v(x). (1) (2) x = [x1, . . . , xn]T is the nonnegative vector of reactant concentrations, and v = [v1, . . . , vm]T is the vector of reaction rates, assumed to be C1. A reaction rate is the rate at which a reaction proceeds to the right and may take any real value. Γ is the (constant) n × m stoichiometric matrix of the reaction system. Since reactant concentrations cannot be negative, it is always reasonable to assume invariance of Rn≥0, i.e. xi = 0 ⇒ x˙ i ≥ 0.
The jth column of Γ , termed Γj , is the reaction vector for the jth reaction, and a stoichiometric matrix is defined only up to an arbitrary signing of its columns. In other words, given any m × m signature matrix D (i.e. any diagonal matrix with diagonal entries ±1), one could replace Γ with Γ D and v(x) with Dv(x). Obviously the dynamical system is left unchanged. The subspace Im(Γ ) of Rn spanned by the reaction vectors is called the stoichiometric subspace. n The intersection of any coset of the Im(Γ ) with R≥0 is called a stoichiometry class.
Two generalisations of (2) which include explicit inflow and outflow of substrates are worth considering. The first of these is a so-called CFSTR x˙ = q(xin − x) + Γ v(x). q ∈ R, the flow rate, is generally assumed to be positive, but we allow q = 0 so that (2) becomes a special case of (3). xin ∈ Rn is a constant nonnegative vector representing the “feed” (i.e., inflow) concentrations. The second class of systems is:
x˙ = xin + Γ v(x) − Q(x).
Here Q(x) = [q1(x1), . . . , qn(xn)]T , with each qi(xi) assumed to be a C1 function satisfying ∂∂xqii > 0, and all other quantities defined as before. Systems (4) include systems (3) with q 6= 0, while systems (2) lie in the closure of systems (4).
Define the m × n matrix V = [Vji] where Vji = ∂∂vxji . A very reasonable,
but weak, assumption about many reaction systems is that reaction rates are
monotonic functions of substrate concentrations as assumed in [
As discussed in [
SR graphs are signed, bipartite multigraphs with two vertex sets VS (termed “Svertices”) and VR (termed “R-vertices”). The edges E form a multiset, consisting of unordered pairs of vertices, one from VS and one from VR. Each edge is signed and labelled either with a positive real number or the formal label ∞. In other words, there are functions sgn : E → {−1, 1}, and lbl : E → (0, ∞) ∪ {∞}. The quintuple (VS , VR, E, sgn, lbl) defines an SR graph.
DSR graphs are similar, but have an additional “orientation function” on their edges, O : E → {−1, 0, 1}. The sextuple (VS , VR, E, sgn, lbl, O) defines a DSR graph. If O(e) = −1 we will say that the edge e has “S-to-R direction”, if O(e) = 1, then e has “R-to-S direction”, and if O(e) = 0, then e is “undirected”. An undirected edge can be regarded as an edge with both S-to-R and R-to-S direction, and indeed, several results below are unchanged if an undirected edge is treated as a pair of antiparallel edges of the same sign. SR graphs can be regarded as the subset of DSR graphs where all edges are undirected.
Both the underlying meanings, and the formal structures, of Petri nets and SR/DSR graphs have some similarity. If we replace each undirected edge in a DSR graph with a pair of antiparallel edges, a DSR graph is simply a Petri net graph, i.e. a bipartite, multidigraph. Similarly, an SR graph is a bipartite multigraph. S-vertices correspond to variables, while R-vertices correspond to processes which govern their interaction. The notions of variable and process are similar to the notions of “place” and “transition” for a Petri net. Edges in SR/DSR graphs tell us which variables participate in each process, with additional qualitative information on the nature of this participation in the form of signs, labels, and directions; edges in Petri nets inform on which objects are changed by a transition, again with additional information in the form of labels (multiplicities) and directions. Thus both Petri net graphs and SR/DSR graphs encode partial information about associated dynamical systems, while neither includes an explicit notion of time.
There are some important differences, however. Where SR/DSR graphs
generally represent the structures of continuous-state, continuous-time dynamical
systems, Petri nets most often correspond to discrete-state, discrete-time
systems, although the translation to a continuous-state and continuous-time
context is possible [
Apart from formal variations between Petri nets and SR/DSR graphs,
differences in the notions of state and time lead naturally to differences in the
questions asked. Most current work using SR/DSR graphs aims to inform on
the existence, nature, and stability of limit sets of the associated local semiflows.
Analogous questions are certainly possible with Petri nets, for example questions
about the existence of stationary probability distributions for stochastic Petri
nets [
SR and DSR graphs can be associated with arbitrary CRNs and more general
dynamical systems [
SR/DSR graphs can be uniquely associated with (2), (3), or (4): in the case of (3) and (4), the inflows and outflows are ignored, and the SR/DSR graph is just that derived from the associated system (2). The construction is most easily visualised via an example. Consider, first, the simple system of two reactions: A + B
C,
A
B This has SR graph, shown in Figure 1, left. If all substrates affect the rates of reactions in which they participate then this is also the DSR graph for the reaction. If, now, the second reaction is irreversible, i.e. one can write A + B
C,
A → B, and consequently the concentration of B does not affect the rate of the second reaction1, then the SR graph remains the same, losing information about irreversibility, but the DSR graph now appears as in Figure 1 right. 1 Note that this is usually, but not always, implied by irreversibility: it is possible for the product of an irreversible reaction to influence a reaction rate. (5) (6)
R1 B 1
A 1 R2 1 C 1 1
R1 B 1 In the usual way, cycles in SR (DSR) graphs are minimal undirected (directed) paths from some vertex to itself. All paths have a sign, defined as the product of signs of edges in the path. Given any subgraph E, its size (or length, if it is a path) |E| is the number of edges in E. Paths of length two will be called short paths. Any path E of even length also has a parity
P (E) = (−1)|E|/2sign(E).
A cycle C is an e-cycle if P (C) = 1, and an o-cycle otherwise. Given a cycle C containing edges e1, e2, . . . , e2r such that ei and e(i mod 2r)+1 are adjacent for each i = 1, . . . , 2r, define: stoich(C) = r Y lbl(e2i−1) − i=1 r Y lbl(e2i) . i=1 Note that this definition is independent of the starting point chosen on the cycle. A cycle with stoich(C) = 0 is termed an s-cycle.
An S-to-R path in an SR graph is a non-self-intersecting path between an
Svertex and an R-vertex. R-to-R paths and S-to-S paths are similarly defined,
though in these cases the initial and terminal vertices may coincide. Any cycle is
both an R-to-R path and an S-to-S path. Two cycles have S-to-R intersection
if each component of their intersection is an S-to-R path. This definition can be
generalised to DSR graphs in a natural way, but to avoid technicalities regarding
cycle orientation, the reader is referred to [
Returning to the family of CRNs in (1), these give SR graphs shown in Figure 2. If all reactants can influence the rates of reactions in which they participate, then these are also their DSR graphs (otherwise some edges may become directed). Each SR graph contains a single cycle, which is an e-cycle (resp. ocycle) if n is odd (resp. even). These cycles all fail to be s-cycles because of the unique edge-label of 2. SYS 1 2 2
SYS 2
SYS 3 2
A function f : X → Rn is injective if for any x, y ∈ X, f (x) = f (y) implies x = y. Injectivity of a vector field on some domain is sufficient to guarantee that there can be no more than one equilibrium on this domain. Define the following easily computable condition on an SR or DSR graph:
Condition (∗): All e-cycles are s-cycles, and no two e-cycles have S-to-R intersection.
Note that if an SR/DSR graph has no e-cycles, then Condition (∗) is trivially
fulfilled. A key result in [
Proof. See Theorem 1 in [
In [
Proof. See Corollary 4.2 in [
Proposition 2 is stronger than Proposition 1 because irreversibility is taken
into account. In the case without outflows (2), attention must be restricted to
some fixed stoichiometric class. The results then state that no stoichiometry
class can contain more than one nondegenerate equilibrium in the interior of the
positive orthant [
A closed, convex, solid, pointed cone K ⊂ Rn is termed a proper cone [
An extremal ray is a one dimensional face of a cone. A proper cone with exactly n extremal rays is termed simplicial. Simplicial cones have the feature that unit vectors on the extremal rays can be chosen as basis vectors for a new coordinate system. Consider some linear subspace A ⊂ Rn. Then any closed, convex, pointed cone K ⊂ A with nonempty interior in A is termed A-proper. If, further, K has exactly dim(A) extremal rays, then K is termed A-simplicial.
Consider some local semiflow φ defined on X ⊂ Rn. Assume that there is some linear subspace A ⊂ Rn with a coset A0 with nonempty intersection with X, and such that φ leaves A0 ∩ X invariant. Suppose further that there is an A-proper cone K such that for all x, y ∈ A0 ∩ X, x > y ⇒ φt(x) > φt(y) for all values of t ≥ 0 such that φt(x) and φt(y) are defined. Then we say that φ|A0 ∩X preserves K, and that φ|A0 ∩X is monotone. If, further, x > y ⇒ φt(x) φt(y) for all values of t > 0 such that φt(x) and φt(y) are defined, then φ|A0 ∩X is strongly monotone. A local semiflow is monotone with respect to the nonnegative orthant if and only if the Jacobian of the vector field has nonnegative off-diagonal elements, in which case the vector field is termed cooperative.
Returning to (3), in the case q = 0, all stoichiometry classes are invariant,
while if q > 0, there is a globally attracting stoichiometry class. Conditions for
monotonicity of φ restricted to invariant subspaces of Rn were discussed
extensively in [
Given a vector y ∈ Rn, define
Q1(y) ≡ {v ∈ Rn | viyi ≥ 0}.
A matrix Γ is R-sorted (resp. S-sorted) if any two distinct columns (resp. rows) Γi and Γj of Γ satisfy Γi ∈ Q1(−Γj ). A matrix Γ 0 is R-sortable (resp. S-sortable) if there exists a signature matrix D such that Γ ≡ Γ 0 D (resp. Γ ≡ DΓ 0 ) is well-defined, and is R-sorted (resp. S-sorted).
Proposition 3. Consider a system of N1C reactions of the form (3) whose stoichiometric matrix Γ is R-sortable, and whose reaction vectors {Γk} are linearly independent. Let S = Im(Γ ). Then there is an S-simplicial cone K preserved by the system restricted to any invariant stoichiometry class, such that each reaction vector is collinear with an extremal ray of K.
Proof. This is a specialisation of Corollary A7 in [
Systems fulfilling the assumptions of Proposition 3, cannot have periodic
orbits intersecting the interior of the positive orthant which are stable on their
stoichiometry class. In fact, mild additional assumptions ensure strong
monotonicity guaranteeing generic convergence of bounded trajectories to equilibria
[
Some more notation is needed for the results to follow. The S-degree (Rdegree) of an SR graph G is the maximum degree of its S-vertices (R-vertices). Analogous to the terminology for matrices, a subgraph E is R-sorted (Ssorted) if each R-to-R (S-to-S) path Ek in E satisfies P (Ek) = 1. Note that E is R-sorted if and only if each R-to-R path Ek of length 2 in E satisfies P (Ek) = 1.
An R-flip on a SR/DSR graph G is an operation which changes the signs on all edges incident on some R-vertex in G. (This is equivalent to exchanging left and right for the chemical reaction associated with the R-vertex). An Rresigning is a sequence of R-flips. An S-flip and S-resigning can be defined similarly. Given a set of R-vertices {Rk} in G, the closed neighbourhood of {Rk} will be denoted G{Rk}, i.e., G{Rk} is the subgraph consisting of {Rk} along with all edges incident on vertices of {Rk}, and all S-vertices adjacent to those in {Rk}.
Proposition 4. Consider a system of N1C reactions of the form (3) with stoichiometric matrix Γ , and whose reaction vectors {Γk} are linearly independent. Define S = Im(Γ ). Associate with the system the SR graph G. Suppose that 1. G has S-degree ≤ 2. 2. All cycles in G are e-cycles.
Then there is an S-simplicial cone K preserved by the system restricted to any invariant stoichiometry class, such that each reaction vector is collinear with an extremal ray of K.
The key idea of the proof is simple: if the system satisfies the conditions of Proposition 4, then the conditions of Proposition 3 are also met. In this case, the extremal vectors of the cone K define a local coordinate system on each stoichiometry class, such that the (restricted) system is cooperative in this coordinate system. This interpretation in terms of recoordinatisation is best illustrated with an example.
Consider SYS 1 from (1) with SR graph shown in Figure 2 left, which can easily be confirmed to satisfy the conditions of Proposition 4. Define the following matrices: Γ , the stoichiometric matrix, has rank 3, and so Proposition 4 applies. Let x1, . . . , x5 be the concentrations of the five substrates involved, v1, v2, v3 be the rates of the three reactions, and vij ≡ ∂∂xvji . Assuming that the system is N1C means that V ≡ [vij ] has sign structure
+ + 0 − 0 sgn(V ) = 0 + + 0 −
− 0 + 0 0 where + denotes a nonnegative quantity, and − denotes a nonpositive quantity. Consider now any coordinates y satisfying x = T y. Note that T is a re-signed version of Γ . Choosing some left inverse for T , say T 0 , gives y1 = x1 − 2x2 + 2x3, y2 = x1 − x2 + 2x3 and y3 = x1 − x2 + x3. (The choice of T 0 is not unique, but this does not affect the argument.) Calculation gives that J = T 0 Γ V T has sign structure
− + + sgn(J ) = + − + , + + − i.e., restricting to any invariant stoichiometry class, the dynamical system for the evolution of the quantities y1, y2, y3 is cooperative. Further, the evolution of {xi} is uniquely determined by the evolution of {yi} via the equation x = T y.
It is time to return to the steps leading to the proof of Proposition 4. In Lemmas 1 and 2 below, G is an SR graph with S-degree ≤ 2. This implies the following: consider R-vertices v, v0 and v00 such that v 6= v0 and v 6= v00 (v0 = v00 is possible). Assume there exist two distinct short paths in G, one from v to v0 and one from v to v00 . These paths must be edge disjoint, for otherwise there must be an S-vertex lying on both A and B, and hence having degree ≥ 3. Lemma 1. Suppose G is a connected SR graph with S-degree ≤ 200, and has some connected, R-sorted, subgraph E containing R-vertices v0 and v . Assume that there is a path C1 of length 4 between v0 and v00 containing an R-vertex not in E. Then either C1 is even or G contains an o-cycle.
Proof. If v0 = v00 , then C1 is not even, then it is itself and e-cycle. Otherwise consider any path C2 connecting v0 and v00 and lying entirely in E. C2 exists since E is connected, and P (C2) = 1 since E is R-sorted. Since G has S-degree ≤ 2, and |C1| = 4, C1 and C2 share only endpoints, v0 and v00 , and hence together they form a cycle C. If P (C1) = −1, then P (C) = P (C2)P (C1) = −1, and so C is an o-cycle. tu Lemma 2. Suppose G is a connected SR graph with S-degree ≤ 2 which does not contain an o-cycle. Then it can be R-sorted.
Proof. The result is trivial if G contains a single R-vertex, as it contains no short R-to-R paths. Suppose the result is true for graphs containing k R-vertices.
Then it must be true for graphs containing k + 1 R-vertices. Suppose G contains k + 1 R-vertices. Enumerate these R-vertices as R1, . . . , Rk+1 in such a way that G− ≡ G{R1,...,Rk} is connected. This is possible since G is connected.
By the induction hypothesis, G− can be R-sorted. Having R-sorted G−, consider Rk+1. If all short paths between Rk+1 and R-vertices in G− have the same parity, then either they are all even and G is R-sorted; or they are all odd, and a single R-flip on Rk+1 R-sorts G. (Note that an R-flip on Rk+1 does not affect the parity of any R-to-R paths in G−.) Otherwise there must be two distinct short paths of opposite sign, between Rk+1 and R-vertices v0 , v00 ∈ G− (v0 = v00 is possible). Since G has S-degree ≤ 2, these paths must be edge-disjoint, and together form an odd path of length 4 from v0 to Rk+1 to v00 . By Lemma 1, G contains an o-cycle. tu
PROOF of Proposition 4. From Lemma 2, if no connected component of G contains an o-cycle then each connected component of G (and hence G itself) can be R-sorted. The fact that G can be R-sorted corresponds to choosing a signing of the stoichiometric matrix Γ such that any two columns Γi and Γj satisfy Γi ∈ Q1(−Γj ). Thus the conditions of Proposition 3 are satisfied. tu 7
Example 1: SYS n from Section 1. It is easy to confirm that the reactions in SYS n have linearly independent reaction vectors for all n . Moreover, as illustrated by Figure 2, the corresponding SR graphs contain a single cycle, which, for odd (even) n is an e-cycle (o-cycle). Thus for even n, Proposition 1 and subsequent remarks apply, ruling out the possibility of more than one positive nondegenerate equilibrium for (2) on each stoichiometry class, or in the case with outflows (4), ruling out multiple equilibria altogether; meanwhile, while for odd n, Proposition 4 can be applied to (2) or (3), implying that restricted to any invariant stoichiometry class the system is monotone, and the restricted dynamical system cannot have an attracting periodic orbit intersecting the interior of the nonnegative orthant.
lowing system of chemical reactions:
A
B,
A
C,
A
D,
B
C
(7)
with SR graph shown in Figure 3. Formally, such systems have R-degree ≤ 2 and
have SR graphs which are S-sorted. Although Proposition 4 cannot be applied,
such “interconversion networks”, with the N1C assumption, in fact give rise to
cooperative dynamical systems [
This example highlights that there is an immediate dual to Lemma 2, and
hence Proposition 4. The following lemma can be regarded as a restatement of
well-known results on systems preserving orthant cones (see [
R1 A
C R3
R4 B
Lemma 3. Let G be an SR graph with R-degree ≤ 2 and containing no o-cycles. Then, via an S-resigning, G can be S-sorted.
Although the S-sorting process is formally similar to the R-sorting one, the interpretation of the result is quite different: changing the sign of the ith row of Γ and the ith column of V is equivalent to a recoordinatisation replacing concentration xi with −xi. Such recoordinatisations give rise to a cooperative system if and only if the original system is monotone with respect to an orthant cone.
sary for monotonicity. Consider the system of three reactions involving four substrates
A with stoichiometric matrix Γ and SR graph shown in Figure 4.
Γ = T = −10 01 00 . Note: i) T has rank 3, ii) Im(Γ ) ⊂ Im(T ), and iii) regarding the columns of T as extremal vectors of a cone K, K has trivial intersection with Im(Γ ). One can proceed to choose some left inverse T 0 of T , and calculate that the Jacobian J = T 0 Γ V T has nonnegative off-diagonal entries. In other words the y-variables define a cooperative dynamical system. The relationship between T and Γ is further discussed in the concluding section.
Note that although K has empty interior in R4, both K and Im(Γ ) lie in the
hyperplane H = Im(T ) defined by x1 + x3 = 0. As K is H-proper, attention can
be restricted to invariant cosets of H. With mild additional assumptions on the
kinetics, the theory in [
tonicity. Consider the following system of 4 chemical reactions on 5 substrates:
A
Γ , the stoichiometric matrix, has rank 3, and the system has SR graph containing
both e- and o-cycles (Figure 5). Further, there are substrates participating in 3
reactions, and reactions involving 3 substrates (and so it is neither R-sortable nor
S-sortable). Thus, all the conditions for the results quoted so far in this paper,
and for theorems in [
R4 A
E
C R1
R3 B
D R2
tonicity. Returning to the system of reactions in (5), the system has SR graph
containing an o-cycle (Figure 1, left). Nevertheless, the system was shown in
[
The results presented here provide only a glimpse of the possibilities for analysis of limit sets of CRNs using graph-theoretic – and more generally combinatorial – approaches. The literature in this area is growing rapidly, and new techniques are constantly being brought into play. Working with the weakest possible kinetic assumptions often gives rise to approaches quite different from those used in the previous study of mass-action systems. Conversely, it is possible that such approaches can be used to provide explicit restrictions on the kinetics for which a system displays some particular behaviour.
The paper highlights an interesting duality between questions of
multistationarity and questions of stable periodic behaviour, a duality already implicit
in discussions of interaction graphs [
One open question regards the relationship between the theory and examples
presented here on monotonicity, and previous results, particularly Theorem 1 in
[
Consider again Examples 3 and 4a. The key fact is that their stoichiometric
matrices admit factorisations Γ = T1T2, taking the particular forms
−1 0 1 1
011 −110 −−110 −100 = −010 001 100 000
0 0 0 −1 0 0 0 1
In each case, the first factor, T1, has exactly one nonzero entry in each row. On
the other hand, the second factor, T2, is S-sorted. The theory in [
A broad open question concerns the extent to which the techniques
presented here extend to systems with discrete-time, and perhaps also
discretestate space. In [
Finally, the work on monotonicity here has an interesting relationship with
examples presented by Kunze and Siegel, for example in [