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:

• John Baez and Jacob Biamonte, Quantum techniques for stochastic physics.

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 CO₂ and HCO₃^-, 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 $S$ of species, a complex of those species is a function $C: S \to \mathbb{N}$.

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

Definition. A reaction network consists of a set $S$ of species, a set $T$ of reactions, and functions $s,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,\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 $\tau \in T$ is labelled by a rate constant $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 $k$ different species. Let $x_i$ be the number of things of the $i$th species. Then the rate equation looks like this:

It’s really a bunch of equations, one for each $1 \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 $i$th species appears as input to this transition $m_i$ times, and as output $n_i$ times. Then the rate equation is

where $r$ 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

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

and an output vector

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

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

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

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:

• John Baez and Jacob Biamonte, Quantum techniques for stochastic physics.

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^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.

$l$ linkage classes. These are indexed by $\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) $\delta_\theta = 0$ , $\theta = 1 , 2 ,..., l$

• (ii) $\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) $\delta_\theta \langle 1, \theta = 1,2 ,..., l$

• (ii) $\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:

• Katalin M. Hangos, Engineering model reduction and entropy-based Lyapunov functions in chemical reaction kinetics, Entropy 12 (2010), 772-797.

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:

$$\frac{d x_i}{d t} = F_i(x) = -x_i f_i(x) + g_i(x) \qquad i = 1, \dots, n$$

where $f_i$ and $g_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:

• Dávid Csercsik, Gábor Szederkényi and Katalin M. Hangos, Parametric uniqueness of deficiency zero reaction networks, Journal of Mathematical Chemistry 50 (2012), 1-8.

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:

• Katalin M. Hangos, Engineering model reduction and entropy-based Lyapunov functions in chemical reaction kinetics, Entropy 12 (2010), 772–797.

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., 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.

with strength $k_1$, $k_2$, respectively.

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

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

• Reactions: There are two reactions. $v_1 \rightarrow v_1'\in R$ and $v_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_k'-v_k$ so it is one dimensional.

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

The general form of the rate equation is given as

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:

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

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

Example: 1-parameter family of equilibrium states

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

The input and output functions then give.

and for $\tau_2$

The general form of the rate equation is

The population of species $A$ ($B$) will be denoted $X_1(t)$ ($X_2(t)$) and the rate of the reaction $\tau_1$ ($\tau_2$) as $r_1$ ($r_2$).

and

The general form of the master equation is

In our case, this becomes

We consider a wave function of the form

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

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

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

The input and output functions then give.

and for $\tau_2$

The general form of the rate equation is

The population of species $A$ ($B$) will be denoted $X_1(t)$ ($X_2(t)$) and the rate of the reaction $\tau_1$ ($\tau_2$) as $r_1$ ($r_2$).

and

The general form of the master equation is

In our case, this becomes

We consider a wave function of the form

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

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

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

Here the rate equation is:

The equilibrium solutions of this equation are given as:

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

The Hamiltonian for the master equation in this example is:

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

By picking $c_1 = X_1$, $c_2 = X_2$ and $c_3 = X_3$. This example has the following Petri net:

Example: a bimolecular reaction

Here the rate equation is:

The Hamiltonian for the master equation is:

From this we arrive at the following

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.

• M. Feinberg (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.

• M. Feinberg (1989). “Necessary and sufficient conditions for detailed balancing in mass action systems of arbitrary complexity”. Chemical Engineering Science 44 (9): 1819–1827. doi:10.1016/0009-2509(89)85124-3.

• M. Feinberg (December 1995). “The existence and uniqueness of steady states for a class of chemical reaction networks”. Archive for Rational Mechanics and Analysis 132 (4): 311–370. doi:10.1007/BF00375614.

• E.D. Sontag and D. Zeilberger. A symbolic computation approach to a problem involving multivariate Poisson distributions. Advances in Applied Mathematics, 44:359-377, 2010.

• Schlosser, P. M. and M. Feinberg (June 1994). “A theory of multiple steady states in isothermal homogeneous CFSTRs with many reactions”. Chemical Engineering Science 49 (11): 1749–1767. doi:10.1016/0009-2509(94)80061-8.

A list of further resources:

• Petri Nets and Ecological Modelling