The Azimuth Project
Chemical reaction network (changes)

Showing changes from revision #34 to #35: Added | Removed | Changed

Idea

Chemical reaction networks are a formalism equivalent to stochastic Petri nets but preferred by chemists. Many deep results on stochastic Petri nets are being proved by chemists—but phrased in terms of chemical reaction networks. They are discussed in detail here:

starting in Chapter 17. This page currently contains a short explanation of the concept of chemical reaction network, or reaction network for short, followed by some quotes from fundamental papers on chemical reaction networks and some examples.

Table of contents

Details

Here’s an example of a reaction network:

This network involves 5 species: that is, different kinds of things. They could be atoms, molecules, ions or whatever: chemists call all of these species, and there’s no need to limit the applications to chemistry: in population biology, they could even be biological species! We’re calling them A, B, C, D, and E, but in applications, we’d call them by specific names like CO2 and HCO3-, or ‘rabbit’ and ‘wolf’, or whatever.

This network also involves 5 reactions, which are shown as arrows. Each reaction turns one bunch of species into another. So, written out more longwindedly, we’ve got these reactions:

A \to B

B \to A

A + C \to D

B + E \to A + C

B + E \to D

In chemistry, a bunch of species is called a ‘complex’. But what is a ‘bunch’, exactly? In a given complex, each species can show up 0,1,2,3… or any natural number of times. So, we can formalize the concept like this:

Definition. Given a set SS of species, a complex of those species is a function C:SC: S \to \mathbb{N}.

In a reaction network, each reaction goes from some complex to another. So, there’s a set TT of reactions, and for each reaction τT\tau \in T a complex s(τ)s(\tau) called its source and a complex t(τ)t(\tau) called its target.

Definition. A reaction network consists of a set SS of species, a set TT of reactions, and functions s,t:T Ss,t: T \to \mathbb{N}^S saying the source and target of each reaction.

One can now check that any reaction network gives a Petri net, and vice versa. In Petri net theory the species are called states or places, while the reactions are called transitions.

In a stochastic Petri net each transition is labelled by a rate constant: that is, a number in (0,)(0,\infty). This lets us write down some differential equations saying how species turn into each other. So, let’s make this definition:

Definition. A stochastic reaction network is a reaction network where each reaction τT\tau \in T is labelled by a rate constant r(τ)(0,)r(\tau) \in (0,\infty).

Every chemical reaction network gives a ‘rate equation’ and a ‘master equation’:

• The master equation says how the probability that we have a given number of things of each species changes with time.

• The rate equation says how the expected number of things of each species changes with time.

The master equation is stochastic: it describes how probabilities change with time. The rate equation is deterministic.

The rate equation

Suppose we have a reaction network with kk different species. Let x ix_i be the number of things of the iith species. Then the rate equation looks like this:

dx idt=??? \frac{d x_i}{d t} = ???

It’s really a bunch of equations, one for each 1ik1 \le i \le k. But what is the right-hand side?

The right-hand side is a sum of terms, one for each transition in our reaction network. So, let’s start by assuming there is just one transition.

Suppose the iith species appears as input to this transition m im_i times, and as output n in_i times. Then the rate equation is

dx idt=r(n im i)x 1 m 1x k m k \frac{d x_i}{d t} = r (n_i - m_i) x_1^{m_1} \cdots x_k^{m_k}

where rr is the rate constant for this transition.

That’s really all there is to it! But we can make it look nicer. Let’s make up a vector

x=(x 1,,x k)[0,) k x = (x_1, \dots , x_k) \in [0,\infty)^k

that says how many things there are of each species. Similarly let’s make up an input vector

m=(m 1,,m k) k m = (m_1, \dots, m_k) \in \mathbb{N}^k

and an output vector

n=(n 1,,n k) k n = (n_1, \dots, n_k) \in \mathbb{N}^k

for our transition. To be cute, let’s also define

x m=x 1 m 1x k m k x^m = x_1^{m_1} \cdots x_k^{m_k}

Then we can write the rate equation for a single transition like this:

dxdt=r(nm)x m \frac{d x}{d t} = r (n-m) x^m

Next let’s do a general reaction network, with lots of transitions. Let’s write TT for the set of transitions and r(τ)r(\tau) for the rate constant of the transition τT\tau \in T. Let n(τ)n(\tau) and m(τ)m(\tau) be the input and output vectors of the transition τ\tau. Then the rate equation is:

dxdt= τTr(τ)(n(τ)m(τ))x m(τ) \frac{d x}{d t} = \sum_{\tau \in T} r(\tau) (n(\tau) - m(\tau)) x^{m(\tau)}

The Deficiency Zero Theorem

The deficiency zero theorem is a major result on the rate equations of chemical networks. A precise statement and proof are given starting in Chapter 17 here:

It was stated in the following paper:

  • Feinberg, M. (1987). Chemical reaction network structure and the stability of complex isothermal reactors: I. The deficiency zero and deficiency one theorems. Chemical Engineering Science 42 (10): 2229–2268. doi:10.1016/0009-2509(87)80099-4.

Here we quote the theorem from this paper.

  • Theorem 6.1.1. (The Deficiency Zero Theorem): For any reaction network of deficiency zero the following statements hold true:

  • (i) If the network is not weakly reversible then, for arbitrary kinetics (not necessarily mass action), the differential equations for the corresponding reaction system cannot admit a positive steady state (i-e. a steady state in P NP^N).

  • (ii) If the network is not weakly reversible then, for arbitrary kinetics (not necessarily mass action), the differential equations for the corresponding reaction system cannot admit a cyclic composition trajectory along which all species concentrations are positive.

  • (iii) If the network is weakly reversible then, for mass action kinetics (but regardless of the positive values the rate constants take), the differential equations for the corresponding reaction system have the following properties: There exists within each positive stoichiometric compatibility class precisely one steady state; that steady state is asymptotically stable; and there is no nontrivial cyclic composition trajectory along which all species concentrations are positive.

The Deficiency One Theorem

In the following, we quote directly from the following paper.

  • Feinberg, M. (1987). Chemical reaction network structure and the stability of complex isothermal reactors: I. The deficiency zero and deficiency one theorems. Chemical Engineering Science 42 (10): 2229–2268. doi:10.1016/0009-2509(87)80099-4.

Recall from Section 2.6 that we can calculate not only the deficiency of a reaction network as a whole but also the deficiency of each of its linkage classes. Recall also that the deficiency of the entire network need not be the same as the sum of the deficiencies computed for the linkage classes separately (Remark 2.6.A).

The deficiency zero networks are precisely those that satisfy both of the following conditions: First, the deficiency of each linkage class is zero; and, second, the deficiency of the entire network is equal to the sum of the deficiencies of the individual linkage classes.

ll linkage classes. These are indexed by θ=1,...,l\theta = 1,...,l

Each linkage class has a deficiency δ θ\delta_\theta

Thus, some of the information given by part (iii) of the Deficiency Zero Theorem can be phrased differently: Consider a mass action system for which the underlying reaction network is weakly reversible and has P linkage classes. Moreover, suppose that the deficiency of the entire network and the deficiencies of the individual linkage classes satisfy the following conditions:

  • (i) δ θ=0\delta_\theta = 0 , θ=1,2,...,l\theta = 1 , 2 ,..., l
  • (ii) θ=1 lδ θ=δ\sum_{\theta = 1}^l \delta_\theta = \delta.

Then, no matter what (positive) values the rate constants take, the corresponding differential equations admit precisely one steady state in each positive stoichiometric compatibility class.

With respect to the preclusion of multiple positive steady states, the Deficiency One Theorem not only relaxes the deficiency requirements given in the statement above, it also replaces the weak reversibility condition by the far milder condition that each linkage class contain no more than one terminal strong linkage class. (Recall that weakly reversible networks are merely very special examples of networks satisfying this milder condition on reaction arrow structure.)

Theorem 6.2.1. (The Deficiency One Theorem): Consider a mass action system for which the underlying reaction network has P linkage classes, each containing just one terminal strong linkage class. Suppose that the deficiency of the network and the deficiencies of the individual linkage classes satisfy the following conditions:

  • (i) δ θ1,θ=1,2,...,l\delta_\theta \langle 1, \theta = 1,2 ,..., l
  • (ii) θ=1 lδ θ=δ\sum_{\theta = 1}^l \delta_\theta = \delta

Then, no matter what (positive) values the rate constants take, the corresponding differential equations can admit no more than one steady state within a positive stoichiometric compatibility class. If the network is weakly reversible, the differential equations admit precisely one steady state in each positive stoichiometric compatibility class.

  • Corollary 6.2.2. A mass action system for which the underlying reaction network has just one linkage class can admit multiple steady states within a positive stoichiometric compatibility class only if the deficiency of the network or the number of its terminal strong linkage classes exceeds one.

The Anderson-Craciun-Kurtz theorem

This paper:

• D. F. Anderson, G. Craciun and T.G. Kurtz, Product-form stationary distributions for deficiency zero chemical reaction networks, arXiv:0803.3042.

proves that for a large class chemical reaction networks, there exist equilibrium solutions of the master equation where the number of things in each state is distributed according to a Poisson distribution. Even more remarkably, these probability distributions are independent, so that knowing how many things are in one state tells you nothing about how many are in another. A nice quote from the paper is:

The surprising aspect of the deficiency zero theorem is that the assumptions of the theorem are completely related to the network of the system whereas the conclusions of the theorem are related to the dynamical properties of the system.

Existence and uniqueness of chemical reaction networks with given rate equations

Given a rate equation, we can ask if there exists a chemical reaction network it comes from, and whether it’s unique. The existence problem has apparently been settled, as discussed here:

She writes:

In the most general case the dynamics of a reaction kinetic system in the concentration (state) space can be described with the following set of ordinary differential equations (ODEs) that originate from the dynamic conservation balance equations constructed for component masses:

dx idt=F i(x)=x if i(x)+g i(x)i=1,,n \frac{d x_i}{d t} = F_i(x) = -x_i f_i(x) + g_i(x) \qquad i = 1, \dots, n

where f if_i and g ig_i are non-negative polynomials, i.e., polynomials with positive coefficients (…) It can be shown that a model (of this type) can always be realized with a reaction kinetic system obeying the mass action law (17).

Here reference (17) is to this book:

  • Hárs, V.; Tóth, J., Qualitative Theory of Differential Equations, chapter On the inverse problem of reaction kinetics, North Holland, 1981, pp. 363–379.

Uniqueness is discussed in this paper:

Abstract: In this paper it is shown that deficiency zero mass action reaction networks containing one terminal linkage class are parametrically and therefore structurally unique with a fixed complex set. Clearly, weakly reversible deficiency zero networks with one linkage class belong to this class. However, it is shown through an illustrative example that deficiency zero networks with several linkage classes can have multiple dynamically equivalent realizations, even if the individual linkage classes are weakly reversible.

The introduction sketches some background:

According to the well-known “fundamental dogma of chemical kinetics” different reaction networks can produce the same kinetic differential equations. This means, that dynamically equivalent representations called realizations - i.e. reaction networks with possibly different structure and/or reaction rate coefficients from the original one - may exist that still lead to the same kinetic differential equations. This fact has a great importance from the viewpoint of analyzing the properties of a reaction kinetic system given by its kinetic differential equations, because some of the most important structural properties, such as (weak) reversibility or deficiency are realization-dependent, i.e. they may change depending of the particular realization.

Reaction kinetic systems form a special sub-class of positive systems with smooth, polynomial nonlinearities in the ordinary differential equation (ODE) description implied by the mass action law. Beside the description of classical chemical reactions, chemical reaction networks (CRNs) are the main building blocks of highly interconnected biochemical networks with complex behavior such as metabolic or cell signalling pathways. Because of the significance of structural properties, it is of great practical and theoretical interest to find realizations of a given reaction kinetic system with desired properties. The first step to this is to solve the so-called inverse problem of reaction kinetics (i.e. the characterization of those polynomial differential equations which are kinetic), that was published in (7). Recently, optimization-based computational algorithms have been presented for the construction of CRN

structures with preferred properties such as reversibility or minimal/maximal number of reactions and complexes in (10) and (11).

Chemical reaction networks and entropy

This paper gives an excellent introduction to the use of a quantity related to entropy in the theory of chemical reaction networks, starting in the subsection ‘the entropy-based Lyapunov function’ in section 3.3:

Robustness

Here are some references on ‘robustness’ for chemical reaction networks:

  • Guy Shinar, Uri Alon, and Martin Feinberg, Sensitivity and robustness in chemical reaction networks,SIAM J. Appl. Math.

    Guy Shinar, Uri Alon, and Martin Feinberg, Sensitivity and robustness in chemical reaction networks, SIAM J. Appl. Math., 69 (2009), 977–998.

    ,69 (2009), 977–998.
  • Guy Shinar and Martin Feinberg, Design principles for robust biochemical reaction networks: What works, what cannot work, and what might almost work, {Mathematical Biosciences](http://www.sciencedirect.com/science/article/pii/S0025556411000307) 231 (2011), 39–48.

Abstract: We bring together recent results that connect the structure of a mass-action reaction network to its capacity for concentration robustness — that is, its capacity to keep the concentration of a critical bio-active species within narrow limits, even against large fluctuations in the overall supply of the network’s constituents.

Examples

Example 5.1, page 13

We consider the following reaction.

S 1S 2 S_1 \rightarrow S_2
S 2S 1 S_2\rightarrow S_1

with strength k 1k_1, k 2k_2, respectively.

  • Chemical Species: There are two chemical species. S={S 1,S 2}S = \{S_1, S_2\}.

  • Complexes: |C|=2|C|=2. v 1v_1 and v 1v_1'. Note the slight abuse of notation, between the species and complexes. It would have been enough to name the complexes S 1S_1 and S 2S_2.

  • Reactions: There are two reactions. v 1v 1Rv_1 \rightarrow v_1'\in R and v 1v 1Rv_1' \rightarrow v_1\in R.

  • Linkage classes: the network is fully connected so there is one linkage class.

  • Reversibility: The network is reversible.

  • Stoichiometric Subspace: Simply given by the span of v kv kv_k'-v_k so it is one dimensional.

  • Deficiency: |C|=2|C|=2. l=1l=1 and s=1s=1 so δ=0\delta = 0.

The general form of the rate equation is given as

ddtx i(t)= τTr(τ)X 1(t) m 1(τ)X k(t) m k(τ)(n i(τ)m i(τ)) \frac{d}{dt}x_i(t) = \sum_{\tau\in T} r(\tau) X_1(t)^{m_1(\tau)}\cdots X_k(t)^{m_k(\tau)}(n_i(\tau)-m_i(\tau))
ddtX 1(t)=k 1X 1+k 2X 2 \frac{d}{dt}X_1(t) = -k_1X_1 + k_2X_2
ddtX 2(t)=k 1X 1k 2X 2 \frac{d}{dt}X_2(t) = k_1X_1-k_2X_2

Examples from Feinberg

This paper:

  • P. M. Schlosser and M. Feinberg, A theory of multiple steady states in isothermal homogeneous CFSTRs with many reactions, Chemical Engineering Science 49 (1994), 1749–1767.

gives some nice examples of chemical reaction networks that have more than one equilibrium state. The simplest one in Table 1 is this:

2A+B3A 2 A + B \to 3 A
3A2A+B 3 A \to 2 A + B

It would be nice to check this one explicitly. This similar-looking example supposedly has just one equilibrium state:

A+2B3A A + 2 B \to 3 A
3AA+2B 3 A \to A + 2 B

By ‘just one equilibrium state’, I presume they mean just a 1-parameter family of equilibrium states. After all, if x A,x Bx_A, x_B is an equilibrium solution of the rate equation, so is cx A,cx Bc x_A, c x_B for any c0c \ge 0.

Example: 1-parameter family of equilibrium states

Here’s a chemical reaction network with a 1-parameter family of equilibrium states:

τ 1:A+2B3A\tau_1: A + 2 B \to 3 A
τ 2:3AA+2B\tau_2: 3 A \to A + 2 B

The input and output functions then give.

m(τ 1)={1,2} m(\tau_1)=\{1,2\}
n(τ 1)={3,0} n(\tau_1)=\{3,0\}

and for τ 2\tau_2

m(τ 2)={3,0} m(\tau_2)=\{3,0\}
n(τ 2)={1,2} n(\tau_2)=\{1,2\}

The general form of the rate equation is

ddtX i= τTr(τ)X m(τ)(n i(τ)m i(τ)) \frac{d}{d t}X_i=\sum_{\tau\in T}r(\tau)X^{m(\tau)}(n_i(\tau)-m_i(\tau))

The population of species AA (BB) will be denoted X 1(t)X_1(t) (X 2(t)X_2(t)) and the rate of the reaction τ 1\tau_1 (τ 2\tau_2) as r 1r_1 (r 2r_2).

ddtX 1=2r 1X 1X 2 22r 2X 1 3 \frac{d}{d t}X_1 = 2 r_1 X_1 X_2^2 - 2 r_2 X_1^3

and

ddtX 2=2r 1X 1X 2 2+2r 2X 1 3 \frac{d}{d t} X_2 = - 2 r_1 X_1 X_2^2 + 2 r_2 X_1^3

The general form of the master equation is

H= τTr(τ)(a n(τ)a m(τ)a m(τ)a m(τ)) H = \sum_{\tau\in T} r(\tau)(a^\dagger^{n(\tau)}a^{m(\tau)} - a^\dagger^{m(\tau)} a^{m(\tau)})

In our case, this becomes

H=r 1(a 1 a 1 )a 1 a 1a 2a 2a 1 a 2 a 2 a 1a 2a 2)+r 2(a 1 a 2 a 2 a 1a 1a 1a 1 a 1 a 1 a 1a 1a 1) H = r_1 (a_1^\dagger a_1^\dagger)a_1^\dagger a_1 a_2 a_2 - a_1^\dagger a_2^\dagger a_2^\dagger a_1 a_2 a_2) + r_2 (a_1^\dagger a_2^\dagger a_2^\dagger a_1 a_1 a_1 - a_1^\dagger a_1^\dagger a_1^\dagger a_1 a_1 a_1 )

We consider a wave function of the form

Ψ:=e bz 1e cz 2 \Psi := e^{b z_1} e^{c z_2}

We then calculate HΨH\Psi. We find solutions where this vanishes, for non-zero Ψ\Psi as

(r 1c 2br 2b 3)z 1 3+(r 2b 3r 1c 2b)z 1z 2 2 (r_1 c^2 b - r_2 b^3)z_1^3 + (r_2 b^3 - r_1 c^2 b)z_1 z_2^2

For this to vanish, each of the terms must vanish separately. Hence we arrive at a system of equations

r 1c 2br 2b 3=0 r_1 c^2 b - r_2 b^3 = 0
r 2b 3r 1c 2b=0 r_2 b^3 - r_1 c^2 b = 0

In this case, and in every case I’ve considered so far, a solution of one, gives a solution of the second.

Example: multiple equilibrium states

2A+B3A 2 A + B \to 3 A
3A2A+B 3 A \to 2 A + B

The input and output functions then give.

m(τ 1)={2,1} m(\tau_1)=\{2,1\}
n(τ 1)={3,0} n(\tau_1)=\{3,0\}

and for τ 2\tau_2

m(τ 2)={3,0} m(\tau_2)=\{3,0\}
n(τ 2)={2,1} n(\tau_2)=\{2,1\}

The general form of the rate equation is

ddtX i= τTr(τ)X m(τ)(n i(τ)m i(τ)) \frac{d}{d t}X_i=\sum_{\tau\in T}r(\tau)X^{m(\tau)}(n_i(\tau)-m_i(\tau))

The population of species AA (BB) will be denoted X 1(t)X_1(t) (X 2(t)X_2(t)) and the rate of the reaction τ 1\tau_1 (τ 2\tau_2) as r 1r_1 (r 2r_2).

ddtX 1=r 1X 1 2X 2r 2X 1 3 \frac{d}{d t}X_1 = r_1 X_1^2 X_2 - r_2 X_1^3

and

ddtX 2=r 1X 1 2X 2+r 2X 1 3 \frac{d}{d t} X_2 = r_1 X_1^2 X_2 + r_2 X_1^3

The general form of the master equation is

H= τTr(τ)(a n(τ)a m(τ)a m(τ)a m(τ)) H = \sum_{\tau\in T} r(\tau)(a^\dagger^{n(\tau)}a^{m(\tau)} - a^\dagger^{m(\tau)} a^{m(\tau)})

In our case, this becomes

H=r 1(a 1 a 1 a 1 a 1a 1a 2a 1 a 1 a 2 a 1a 1a 2)+r 2(a 1 a 1 a 2 a 1a 1a 1a 1 a 1 a 1 a 1a 1a 1) H = r_1 (a_1^\dagger a_1^\dagger a_1^\dagger a_1 a_1 a_2 - a_1^\dagger a_1^\dagger a_2^\dagger a_1 a_1 a_2) + r_2 (a_1^\dagger a_1^\dagger a_2^\dagger a_1 a_1 a_1 - a_1^\dagger a_1^\dagger a_1^\dagger a_1 a_1 a_1 )

We consider a wave function of the form

Ψ:=e bz 1e cz 2 \Psi := e^{b z_1} e^{c z_2}

We then calculate HΨH\Psi. We find solutions where this vanishes, for non-zero Ψ\Psi as

(r 1b 2cr 2b 3)z 1 3+(r 2b 3r 1b 2c)z 1 2z 2 (r_1 b^2 c - r_2 b^3)z_1^3 + (r_2 b^3 - r_1 b^2 c)z_1^2 z_2

For this to vanish, each of the terms must vanish separately. Hence we arrive at a system of equations

r 1b 2cr 2b 3=0 r_1 b^2 c - r_2 b^3 = 0
r 2b 3r 1b 2c=0 r_2 b^3 - r_1 b^2 c = 0

In this case a solution of one of these equations gives a solution of the other.

A three concentration example

We consider now a more complicated example from Feinberg (page 2253). From equation 6.7 we have that

r 1:A 12A 1 r_1: A_1 \rightarrow 2 A_1
r 2:2A 1A 1 r_2: 2 A_1 \rightarrow A_1
r 3:A 1+A 2A 3 r_3: A_1 + A_2 \rightarrow A_3
r 4:A 3A 1+A 2 r_4: A_3 \rightarrow A_1 + A_2
r 5:A 32A 2 r_5: A_3 \rightarrow 2 A_2
r 6:2A 2A 3 r_6: 2 A_2 \rightarrow A_3

Here the rate equation is:

ddtX 1=r 1X 1r 3X 1X 2r 2X 1 2+r 4X 3 \frac{d}{d t} X_1 = r_1 X_1 - r_3 X_1 X_2 - r_2 X_1^2 + r_4 X_3
ddtX 2=r 3X 1X 2+(r 4+2r 5)X 32r 6X 2 2 \frac{d}{d t} X_2 = -r_3 X_1 X_2 + (r_4+2 r_5) X_3 - 2 r_6 X_2^2
ddtX 3=r 3X 1X 2+r 6X 2 2(r 4+r 5)X 3 \frac{d}{d t} X_3 = r_3 X_1 X_2 + r_6 X_2^2 - (r_4+r_5) X_3

The equilibrium solutions of this equation are given as:

X 1=r 1r 2 X_1 = \frac{r_1}{r_2}
X 2=r 1r 3r 5r 2r 4r 6 X_2 = \frac{r_1 r_3 r_5}{r_2 r_4 r_6}
X 3=r 1 2r 3 2r 5r 2 2r 4 2r 6 X_3 = \frac{r_1^2 r_3^2 r_5}{r_2^2 r_4^2 r_6}

These concentration levels make the rate equations zero for all time.

The Hamiltonian for the master equation in this example is:

H=r 1(a 1 a 1 a 1a 1 a 1)+r 2(a 1 a 1a 1a 1 a 1 a 1a 1)+r 3(a 3 a 1a 2a 1 a 2 a 1a 2)+r 4(a 1 a 2 a 3a 3 a 3)+r 5(a 2 a 2 a 3a 3 a 3)+r 6(a 3 a 2a 2a 2 a 2 a 2a 2) H = r_1(a_1^\dagger a_1^\dagger a_1 - a_1^\dagger a_1) + r_2(a_1^\dagger a_1 a_1 - a_1^\dagger a_1^\dagger a_1 a_1)+ r_3(a_3^\dagger a_1 a_2 - a_1^\dagger a_2^\dagger a_1 a_2) + r_4(a_1^\dagger a_2^\dagger a_3 - a_3^\dagger a_3) + r_5(a_2^\dagger a_2^\dagger a_3 - a_3^\dagger a_3) + r_6(a_3^\dagger a_2 a_2 - a_2^\dagger a_2^\dagger a_2 a_2)

We let this Hamiltonian act on Ψ:=e c 1z 1e c 2z 2e c 3z 3\Psi := e^{c_1 z_1}e^{c_2 z_2}e^{c_3 z_3} and to satisfy HΨ=0H\Psi = 0 there must exist a choice of c 1c_1, c 2c_2 and c 3c_3 that will cause all the following equations to vanish together.

z 1;0=r 1c 1+r 2c 1 2 z_1; 0 = -r_1c_1 + r_2 c_1^2
z 1 2;0=r 1c 1r 2c 1 2 z_1^2; 0 = r_1 c_1 - r_2 c_1^2
z 3;0=r 3c 1c 2(r 4+r 5)c 3+r 6c 2 2 z_3; 0 = r_3 c_1 c_2 - (r_4 +r_5) c_3 + r_6 c_2^2
z 2 2;0=r 6c 2 2+r 5c 3 z_2^2; 0 = -r_6 c_2^2 + r_5 c_3
z 1z 2;0=r 4c 3r 3c 1c 2 z_1 z_2; 0 = r_4 c_3 - r_3 c_1 c_2

By picking c 1=X 1c_1 = X_1, c 2=X 2c_2 = X_2 and c 3=X 3c_3 = X_3. This example has the following Petri net:

Example: a bimolecular reaction

k 1:2X 1+X 22X 3 k_1 : 2 X_1 + X_2 \rightarrow 2 X_3
k 2:2X 32X 1+X 2 k_2: 2 X_3 \rightarrow 2 X_1 + X_2

Here the rate equation is:

ddtX 1=2k 1X 1 2X 2+2k 2X 3 2 \frac{d}{d t} X_1 = -2 k_1 X_1^2 X_2 + 2 k_2 X_3^2
ddtX 2=k 1X 1 2X 2+k 2X 3 2 \frac{d}{d t} X_2 = - k_1 X_1^2 X_2 + k_2 X_3^2
ddtX 3=2k 1X 1 2X 2+2k 2X 3 2 \frac{d}{d t} X_3 = 2 k_1 X_1^2 X_2 + -2 k_2 X_3^2

The Hamiltonian for the master equation is:

H=k 1(a 3 a 3 a 1a 1a 2a 1 a 1 a 1a 1a 2 a 2)+k 2(a 1 a 1 a 2 a 3 a 3 a 3a 3) H = k_1(a_3^\dagger a_3^\dagger a_1 a_1 a_2 - a_1^\dagger a_1^\dagger a_1 a_1 a_2^\dagger a_2) + k_2 (a_1^\dagger a_1^\dagger a_2^\dagger -a_3^\dagger a_3^\dagger a_3 a_3)

From this we arrive at the following

0=k 1c 1 2c 2k 2c 3 2 0 = k_1 c_1^2 c_2 - k_2 c_3^2
0=k 2c 3 2k 1c 1 2c 2 0 = k_2 c_3^2 - k_1 c_1^2 c_2

References

Here we list some key references on chemical reaction networks. Additional references can be found in other pages and in the blog articles.

  • D. F. Anderson, G. Craciun and T.G. Kurtz,

Product-form stationary distributions for deficiency zero chemical reaction networks, arXiv:0803.3042.

A list of further resources: