Chemical Reaction Network Theory

Chemical Reaction Network Theory

Chemical reaction network theory (CRNTheory) developed by Feinberg, Horn and Jackson [1,2], connects qualitative properties of ordinary differential equations corresponding to a reaction network to the network structure. It gives information on whether or not multistationarity is possible and, if possible, then using CRN Toolbox a pair of steady states and the corresponding parameter vector i.e. reaction rate constants can be obtained. In particular, its assertions are independent of specific parameter values.The only assumption is that all kinetics is of mass-action form. The theory introduces new concepts, such as deficiency, linkage classesetc. of a reaction network, and gives conditions on such networks for the existence, uniqueness, multiplicity (multiequilibria, multiple steady states) and stability of steady states (equilibrium points). What follows is adapted from Feinberg [1].

Reaction Networks and Their Structure

Complexes of a Network. The complexes of a network are the entities that appear before and after the reaction arrows. For example, in the following chemical reaction network:

The complexes are {A1 + A2, A3,A4 + A5, 2A6,A1 + A6, A7}, and the number n of the complexes of this network is six.

The Reaction Vectors for a Network. Associated with each reaction of an N species network there is a reaction vector in N dimensional Euclidean space IRN formed by subtracting the "reactant" complex vector from the "product" complex vector. For example, the full set of six reaction vectors for the network above is

{ e3 - e1 - e2, e4 + e5 - e3,2e6 - e3, e3 - 2e6, e7 - e1 - e6, e1 + e6 - e7} (2.1)

where ei is a unit base of IRNcorresponding to ith species and has one in the ith place and zero in all other places.

The Stoichiometric Subspace for a Network. The stoichiometric subspace of an N species network is the linear subspace comprised of all possible linear combinations of the reaction vectors set. The dimension s of the stoichiometric subspace is the number of elements in the largest linearly independent set composed of reaction vectors for the network. For the network given above the reaction vectors { e3 - e1 - e2, e4 + e5 - e3,2e6 - e3, e1 + e6 - e7} are linearly independent. Therefore, s = 4. i.e. among six vectors (reactions), four vectors (reactions) are linearly independent.

The Linkage Classes of a Network. A linkage class is a set of complexes which are all linked to one another either directly or indirectly such that no complex in the setis linked to any outside the set. The number of linkage classes lfor the example mechanism aboveis two, since there are two disconnected pieces in the standard reaction diagram.

Weakly Reversible Networks. The network is weakly reversible, whenever there is a directed arrow pathway (consisting of one or more arrows) pointing from one complex to another, and there is a directed arrow pathway from the second complex to the first. In the following network since there is a directed arrow path connecting 2A1 back to A3 (via A2), this network is weakly reversible.

.

The Deficiency of a Network. Deficiency is an non-negative integer δ, and is defined as, where n represents the number of complexes of a network, l the number of its linkage classes, and s the dimension of its stoichiometric subspace.

The Deficiency Zero Theorem. If the deficiency of a reaction mechanism is zero,then the following statements are true:

• If the mechanism is not weakly reversible, then for arbitrary kinetics (mass action or otherwise), the induced dynamic equations cannot give rise to an equilibrium in V+ where V+ is the set of all vectors of IRN which have no negative components (since concentrations are nonnegative entities).
• If the mechanism is not weakly reversible, then for arbitrary kinetics (mass action or otherwise), the induced dynamic equations cannot give rise to a composition cycle in V+ containing a point in V+.
• If the mechanism is weakly reversible, then for the mass action kinetics with any choice of positive rate constants there exists in each stoichiometric compatibility class of V+(vectorswhich are the intersections of a parallel of stoichiometric subspace with V+)exactly one equilibrium composition, every equilibrium composition in V+ is asymptotically stable. Relative to the stoichiometric compatibility class in which it resides, and the dynamic equations cannot give rise to a composition cycle in V+. In short, weakly reversible mechanism with zero deficiency gives rise to monostability i.e. one equilibrium point.

The Deficiency One Theorem. Consider a weakly reversible network with l linkage classes. Let δ denotes the deficiency of the network; let δθ denotes the deficiency of the δθth linkage class, θ=1,2,..,l; and suppose that both of the following conditions are satisfied:

and

For a weakly reversible network,if there is only one linkage class, then the induced differential equations can admit multiple equilibria within a positive stoichiometric class only if, the deficiency of the network is two or more.

The Advanced Deficiency Theory. The term advanced deficiency is used when deficiency of a network is greater than one. The advanced deficiency theory has the capability to analyze any network that falls outside the scope of deficiency one theory. As with deficiency one theory, the advanced deficiency theory will determine if a network does not have the qualitative capacity to support multiple steady states. Addition of the advanced deficiency theory one can now use reaction network theory to analyze a greater range of mechanisms, including multiple-pathway mechanisms.

The above summarized CRN theory can readily be applied to a given reaction mechanism using chemical reaction network toolbox (though restricted to twenty complexes) that is developed by Feinberg (

Model Equations

Model 1(mitochondria-dependent apoptosis model in the absence of NO effects): The biochemical reactions and their corresponding ODEs obtained from the assumed mass-action kinetics are given in Tables 1 and 2. Table 3 gives the reaction rate constants; Table 4 gives the formation and degradation rate equations and Table 5 gives the parameter values;Table 6 gives the initial species concentrations; Tables 7 and 8 give the initial concentrations, in the caspase-9 activation module, to obtain steady states I and II.

Model 2 (NO effects on the apoptotic pathway): Table 9 gives NO pathway reactions; Table 10 gives ODEs obtained from assumed mass-action kinetics; Table 11 gives reaction rate constants.

Model 3 (overall model combining models 1 and 2): Table 12 gives the bridging reactions; Table 13 gives the bridging equations used to form model 3; Table 14 gives the rate constants of the bridging reactions.

In these Tables the term J denotes the net formation rate of a species due to a reaction (defined as the forward minus reverse reaction rate).For example, in Table 2 for equation index (1), J1b means the following: In Table 1, the net formation rate for the reaction index (1b) is where the subscript in r1bp denotes forward (plus) reaction rate and r1bm denotes reverse (minus) reaction rate of reaction (1b).

Note that to save space c3z, c3a (or casp3), c8a (or casp8), c9z, c9a (or casp9) abbreviations are used to represent procaspase-3, caspase-3, caspase-8, procaspase-9 and caspase-9 respectively in the following tables.

Table 1. Reactions of the mitochondria-dependent apoptosis model (model 1)

Reaction Index / Reactions of Model I / Forward Reaction
Rate Expressions / Forward Reaction
Rate Expressions
(1) / / r1p=k1p[Cc][Ap] / r1m=k1m[CcAp]
(1b) / / r1bp=k1bp[CcAp]p / r1bm=k1bm[Apop]
(2) / / r2p=k2p[Apop][c9z] / r2m=k2m[Apopc9z]
(3) / / r3p=k3p[Apopc9z][c9z] / r3m=k3m[Apopc9z2]
(3f) / / r3f=k3f[Apopc9z2]
(4) / / r4p=k4p[Apopc9a2] / r4m=k4m[Apopc9a][c9a]
(4b) / / r4bp=k4bp[Apopc9a] / r4bm=k4bm[Apop][c9a]
(15) / / r15p=k15p[Apop][HSP70] / r15m=k15m[ApopHSP]
(5) / / r5p=k5p[c9a][IAP] / r5m=k5m[IAP9]
(5b) / / r5bp=k5bp[Apopc9a][IAP] / r5bm=k5bm[IAPA9]
(5c) / / r5cp=k5cp[Apopc9a2][IAP] / r5cm=k5cm[IAPA29]
(6) / / r6p=k6p[c3z][c9a] / r6m=k6m[c93]
(6f) / / r6f=k6f[c93]
(6b) / / r6bp=k6p[c3z][Apopc9a2] / r6bm=k6m[cA93]
(6bf) / / r6bf=k6f[cA93]
(7) / / r7p=k7p[c3a][IAP] / r7m=k7m[IAP3]
(8) / / r8p=k8p[c8a][Bid] / r8m=k8m[c8B]
(8f) / / r8f=k8f[c8B]
(8) / / r8pp=k8p[c3a][Bid] / r8mp=k8m[c3B]
(8f) / / r8fp=k8f[c3B]
(9) / / r9p=k9p[c3a][Bcl2] / r9m=k9m[c3L]
(9f) / / r9f=k9f[c3L]
(11) / / r11=k11[tBid]
(12a) / / r12a=k12a[tBidmito][Bax]
(12b) / / r12b=k12b[tBidBax][Bax]
(13) / / r13=k13[Bcl2][Bax]
(14) / / r14=k14[Bax2][Ccmito]

Table 2. Mathematical expression of model 1 in terms of ordinary differential equations

Equation Index / Ordinary Differential Equations (Model 1)
(1) / d[Apop]/dt = J1b - J2 + J4b - J15 + jApop
(2) / d[c9z]/dt = -J2 - J3 + Jp9
(3) / d[Apopc9z]/dt = J2 - J3
(4) / d[Apopc9z2]/dt = J3 - J3f
(5) / d[Apopc9a2]/dt = J3f - J4 - J5c - J6b + J6bf
(6) / d[Apopc9a]/dt = J4 - J4b - J5b
(7) / d[c9a]/dt = J4 + J4b - J5 - J6 + J6f + jp9a
(8) / d[c3a]/dt = J6f + J6bf - J7 - J8 + J8f - J9 + J9f - μ×c3a
(9) / d[HSP70]/dt = -J15 + jHsp
(10) / d[ApopHSP]/dt = J15 + jAph
(11) / d[Cc]/dt = j14 - J1 - μ×Cc + MPTP×Ccmito
(12) / d[Ap]/dt = -J1 + JAp
(13) / d[CcAp]/dt = J1 - 7×J1b
(14) / d[IAP]/dt = -J5 - J5b - J5c - J7 + JIAP
(15) / d[IAP9]/dt = J5
(16) / d[IAPA9]/dt = J5b
(17) / d[IAPA29]/dt = J5c
(18) / d[c3z]/dt = -J6 - J6b + Jp3 - 10×nop×c3z
(19) / d[c93]/dt = J6 - J6f
(20) / d[cA93]/dt = J6b - J6bf
(21) / d[IAP3]/dt = J7
(22) / d[c8a]/dt = Jp8 - J8 + J8f
(23) / d[Bid]/dt = -J8 -J8p + jbidp
(24) / d[c8B]/dt = J8 - J8f
(25) / d[tBid]/dt = J8f + J8fp - j11 + j12b -μ×tBid + tbid0
(26) / d[c3B]/dt = J8p - J8fp
(27) / d[Bcl2]/dt = -J9 + jbcl2p - j13
(28) / d[c3L]/dt = J9 - J9f
(29) / d[tBidmito]/dt = j11 - j12a - μ×tBidmito
(30) / d[tBidBax]/dt = j12a - j12b -μ×tBidBax
(31) / d[Bax]/dt = jbax - j12a - j12b - j13
(32) / d[Bax2]/dt = j12b -μ×Bax2
(33) / d[Ccmito]/dt = jCcmito - j14 - MPTP×Ccmito

Table 3. Reaction rate constants of the mitochondria-dependent apoptosis model (model 1)

Reaction
Index / Reaction of Model I / Forward Reaction
Rate Constants / Backward Reaction
Rate Constants
(1) / / /
(1b) / / /
(2) / / /
(3) / / /
(3f) / /
(4) / / /
(4b) / / /
(15) / / /
(5) / / /
(5b) / / /
(5c) / / /
(6) / / /
(6f) / /
(6b) / / /
(6bf) / /
(7) / / /
(8) / / /
(8f) / /
(8) / / /
(8f) / /
(9) / / /
(9f) / /
(11) / /
(12a) / /
(12b) / /
(13) / /
(14) / /

Table 4. Formation and degradation rate equations of model 1

JAp = 0.0001×a1 - μ×Ap
JIAP = 0.0001×a2 - μ×IAP
Jp3 = 0.1×0.0001×a3 - μ×c3z
Jp8 = -μ×c8a
jbidp = 0.0001×a5 - μ×Bid
jbcl2p = 0.0001×a6×p53thresh4/(p534 + p53thresh4) - μ×Bcl2
jbax = 0.0001×a7×(1 + p534/(p534 + p53thresh4)) - μ×Bax
jccmito = 0.0001×a8 - μ×Ccmito
Jp9 = 0.00054366
jAph = -0.0001×ApopHSP
jHsp = 0.0001
jp9a= -0.0001×c9a
jApop = 0.0001

Table 5. Parameters used in simulation of model 1

Parameters In Model I
(1) / a1 = 3
(2) / a2 = 0.3
(3) / a3 = 3
(4) / a5 = 0.3
(5) / a6 = 0.8
(6) / a7 = 0.3
(7) / a8 = 3
(8) / a9 = 1
(9) / a10 = 1
(10) / a11 = 1
(11) / a12 = 1
(12) / a13 = 3
(13) / μ = 0.002×a13
(14) / P53thresh = 0.004
(15) / cooperativity (p) = 1
(16) / nop = 0
(17) / tBid0 = 0
(18) / p53 = 0.0066
(19) / MPTP = 0

Table 6. Initial concentrations of components of model 1

Initial Concentrations
[CcAp]0 = 0 M
[Ap]0 = 0.004 M
[Cc]0 = 0 M
[IAP]0 = 0.004 M
[IAP9]0 = 0 M
[IAPA9]0 = 0 M
[IAPA29]0 = 0 M
[c9a]0 = 0 M
[c3z]0 = 0.004 M
[c3a]0 = 0.00001 M
[cA93]0 = 0 M
[IAP3]0 = 0 M
[Bid]0 = 0.004 M
[c8a]0 = 0.00001 M
[c8B]0 = 0 M
[c3B]0 = 0 M
[Bcl2]0 = 0.004 M
[c3L]0 = 0 M
[Bax]0 = 0.004 M
[Ccmito]0 = 0.004 M

Table 7. Initial concentrations for obtaining steady state I

[Apop]0 = 0.15 M
[c9z]0 = 2.32 M
[Apopc9z]0 = 0.31 M
[Apopc9z2]0 = 1 M
[Apopc9a2]0 = 0.58M
[Apopc9a]0 = 0.31M
[c9a]0 = 5.43M
[HSP70]0 = 3.15 M
[ApopHSP]0 = 1M

Table 8. Initial concentrations for obtaining steady state II

[Apop]0 = 3.16M
[c9z]0 = 0.31M
[Apopc9z]0 = 2.32M
[Apopc9z2]0 = 1M
[Apopc9a2]0 = 1.59M
[Apopc9a]0 = 2.33M
[c9a]0 = 5.43M
[HSP70]0 = 0.15 M
[ApopHSP]0 = 1M

Table 9. NO pathway reactions affecting on apoptotic pathway (model 2)

Reaction Index / Reactions of Model 2 / Reaction Rate Expressions
(1) / /
(2) / /
(3) / /
(4) / /
(5) / /
(6) / /
(7) / /
(8) / /
(9) / /
(10) / /
(11) / /
(12a) / /
(13) / /
(m) / /
(12bp) (12bm) / /

(14) / /
(15) / /
(16) / /
(17) / /
(17b) / /

Table 10. Mathematical expression of model 2 in terms of ordinary differential equations

Equation Index / Ordinary Differential Equations of Model 2
(1) / d[NO]/dt = r1NO – r4NO – 2r12aNO – r12bpNO + r12bmNO + 2r14NO – r15NO – r16NO
(2) / d[O2-]/dt = r2NO – r4NO – r5NO – r10NO – 2r17Bno
(3) / d[ONOO-]/dt = r4NO – r6NO – r7NO – r8NO – r9NO
(4) / d[GSH]/dt = r3NO – r6NO – r11NO + 2rm – r17NO – 2r17bNO
(5) / d[GSNO]/dt = r6NO – 2r10NO + r11NO – 2r14NO + r17NO
(6) / d[N2O3]/dt = – r11NO + r12bpNO – r12bmNO – r13NO
(7) / d[NO2]/dt = 2r12aNO – r12bpNO + r12bmNO
(8) / d[CcOx]/dt = – r15NO
(9) / d[FeLn]/dt= – r16NO + r17NO

Table 11. Reaction rate constants of model 2

Reaction Index / Reactions of Model 2 / Reaction Rate Constants
(1) / /
(2) / /
(3) / /
(4) / /
(5) / /
(6) / /
(7) / /
(8) / /
(9) / /
(10) / /
(11) / /
(12a) / /
(13) / /
(m) / /
(12bp) (12bm) / /

(14) / /
(15) / /
(16) / /
(17) / /
(17b) / /

Table 12. Additional reactions combining model 1 and model 2 (model 3)

Reaction Index / Reactions of Combining Model I and Model 2 / Reaction Rate Constants
(1) / / r18NO = k18NO[ONOO-][MPTPc1]
(2) / / r19NO = k19NO[N2O3][casp8]
(3) / / r20NO = k20NO[FeLnNO][casp8]
(4) / / r21NO = k21NO[FeLnNO][casp9]
(5) / / r22NO = k22NO[FeLnNO][casp3]

Table 13. Bridging equations used to form model 3

Equation Index / Ordinary Differential Equations of Model 3
(1) / d[ONOO-]/dt = r4NO – r6NO – r7NO – r8NO – r9NO – r18NO
(2) / d[MPTPc1]/dt = – r18NO
(3) / d[N2O3]/dt = – r11NO + r12bpNO – r12bmNO – r13NO – r19NO
(4) / d[c8a]/dt = Jp8 - J8 + J8f – r19NO – r20NO
(5) / d[FeLnNO]/dt = r16NO – r17NO – r20NO – r21NO – r22NO
(6) / d[FeLn]/dt = – r16NO + r17NO + r20NO + r21NO + r22NO
(7) / d[c9a]/dt = J4 + J4b - J5 - J6 + J6f + jp9a – r21NO
(8) / d[c3a]/dt = J6f + J6bf - J7 - J8 + J8f - J9 + J9f - μ×c3a – r22NO
(9) / d[Cc]/dt = j14 – J1 – μ×Cc + k[MPTP][Ccmito] where k = 1 M-1s-1

Table 14. Reaction rate constants forthe bridging reactions when model 3is formed

Reaction Index / Reactions of Combining Model I and Model 2 / Reaction Rate Constants
(1) / /
(2) / /
(3) / /
(4) / /
(5) / /

Methodology

The bistable caspase-9 activationmodulein model 1 found by the CRN Toolbox is given in Table 16. The normalized reaction rate constants are given in Table 17 and the normalized steady state concentrations are given in Table 18. Tables 19 and 20 give the steady states obtained when the reaction rate constants as obtained from CRN Toolbox are used in XPPAUT. The difference is due to round off errors. For example for apop (denoted by A) the CRN Toolbox output is 0.15718708and 0.1568525 for the XPPAUT output for steady state I.

Table15. Representation of the species and complexes in the caspase-9 activation module

apop  A / apop.(pro9)2 D / casp9 G
pro9  B / apop.(casp9)2 E / HSP70  H
apop.pro9  C / apop.casp9  F / apop. HSP70  I

Table16. Representation of the reactions in the caspase-9 activation module

Table 17. Rate constants of the reaction network of the caspase-9 activation module

Table 18. Steady state concentrations for the caspase-9 activation module reaction network

Table19. Steady state I of caspase-9 activation module

Table20. Steady state II of caspase-9 activation module

Dimensionalizing the Caspase-9 Activation Module

The reaction rate constants given by CRNToolbox are normalized reaction rate constants. Therefore, order of magnitude of molarity unit of concentrations can be assigned arbitrarily to these rate constants.The reaction rate constants obtained from CRN Toolbox are also checked, herein, for the caspase-9 activation module to see if they are biologically meaningful i.e. whether they have the same order of magnitudes with those available in the literature. Table 21 lists reaction rate constants given by CRN Toolbox, Bagci et al., and the corresponding dimensionalized rate constants.

Table 21. Comparison of reaction rate constants given by CRN Toolbox Bagci et al. (2008) and dimensionalized reaction rate constants used for caspase-9 activation module

Reaction / Rate Constant by CRNT
(dimensionless) / Rate Constant by Bagci et al. (2008)
(in terms of μM) / Reaction Rate Constant
(in terms of μM)
/ Apf=1 / - / Apf=1
/ pro9f=5.44 / pro9f=0.0003 / pro9f=0.0003
/ HSP70f=1 / - / HSP70f=1
/ k15p=4.03
k15m=1.00 / - / k15p=4.03
k15m=1.00
/ k2p=10.23
k2m=3.19 / k2p=10
k2m=0.5 / k2p=10.23
k2m=3.19
/ k3p=5.14
k3m=1.00 / k3p=10
k3m=0.5 / k3p=5.14
k3m=1.00
/ k3f=2.72 / k3f=0.1 / k3f=2.72
/ k4p=6.39
k4m=0.59 / k4p=5
k4m=0.5 / k4p=6.39
k4m=0.59
/ k4bp=13.05
k4bm=1.60 / k4bp=5
k4bm=0.5 / k4bp=13.05
k4bm=1.60
/ apop. HSP70deg=1 / - / apop. HSP70deg=1
/ p9deg=1 / p9deg=0.006 / p9deg=0.006

SR Graphs

In what follows are explained very briefly how a reaction network gives rise to species-reaction (SR) graph and are taken almost verbatim from Craciun et al. [3].Further details can be found this paper. The reactions taken place in a reaction network are written in reaction symbols which are rectangles and each species are written in species symbols which are circles. If a species appears within a reaction, then an arc is drawn from the species symbol to the reaction symbol and the arc is labeled with the name of the complex in which the species appears. The following terminology is needed to infer bistability.

Complex pair (c-pair) in the SR graph: This means a pair of arcs that are adjacent to the same reaction symbol and that carry the same complex label.

A cycle: This means a closed path in which no arc or vertex is traversed twice. A cycle can be an odd-cycle, even-cycle and 1-cycle and these classifications are not mutually exclusive; for example an odd-cycle can also be a 1-cycle. A cycle containing an odd number of c-pairs is called an odd-cycle, and a cycle containing even number of c-pairs is called an even-cycle. If a cycle is such that the stoichiometric coefficient associated with each of its arcs is a 1. Two cycles split a c-pair if both arcs of the c-pair are among the arcs of the two cycles and one of the arcs is contained in just one of the cycles.

Theorem: Consider a reaction network for which the SR graph has both of the following properties.

i)Each cycle in the graph is a 1-cycle, an odd-cycle, or both

ii) No c-pair is split by two even-cycles. For such reaction network, the corresponding mass-action differential equations cannot admit more than one positive steady state, no matter positive what values the rate constants, effluent coefficients, or species supply rates take.

When the conditions i and ii are not satisfied the theorem gives no information. To say that a network does have the capacity for bistability means that positive rate constants were found to elicit bistability.

Using the notation given in Table 16 the SR graph of the caspase-9 module is given in the Figure below. This module is already found to be bistable using CRN Toolbox. The cycle no. 1 consisting of the species ACDEFA in the SR graph is an even cycle because it does not have any c-pairs. The cycle no.2 consisting the species ABDEGA is also an even cycle because it has two c-pairs e.g., A+B and A+G. These two cycles split a c-pair, namely F+G because cycle no.1 includes one arc of F+G c-pair and the cycle no.2 includes the other arc of the c-pair. For this SR graph, the condition ii is violated. Therefore, to elicit the bistability of this reaction network, positive parameters (rate constants) have to be found.

REFERENCES

[1] M. Feinberg (1977) in N. Amundson and L. Lapidus (Ed.), Mathematical Aspects of Mass-Action Kinetics, Ch.1 in Chemical Reaction Theory: A Review, Prentice-Hall, Englewood Cliffs.

[2] F. Horn and R. Jackson (1972), General Mass Action Kinetics, Archive Rational Mech. Anal. 47, 81-116

[3] G. Craciun, Y. Tang and M. Feinberg, (2006), Understanding Bistability in Complex Enzyme-Driven Reaction Networks, Proc. Natl. Acad. Sci. USA, 103 (23), 8697-8702.

1