The Azimuth Project
Blog - network theory (part 18)

This page is a blog article in progress, written by John Baez and Jacob Biamonte. To see discussions of this article as it was being written, go to the Azimuth Forum. For the final polished version, visit the Azimuth Blog.

joint with Jacob Biamonte

Last time we explained how ‘reaction networks’, as used in chemistry, are just another way of talking about Petri nets. We stated an amazing result on reaction networks: Feinberg’s deficiency zero theorem. This settles quite a number of questions about chemical reactions. Now let’s illustrate it with an example.

Our example won’t show how powerful this theorem is: it’s too simple. But it’ll help explain the ideas involved.

Diatomic molecules

A diatomic molecule consists of two atoms of the same kind, stuck together:

At room temperature there are 5 elements that are diatomic gases: hydrogen, nitrogen, oxygen, fluorine, chlorine. Bromine is a diatomic liquid, but easily evaporates into a diatomic gas:

Iodine is a crystal at room temperatures:

but if you heat it a bit, it becomes a diatomic liquid and then a gas:

so people often list it as a seventh member of the diatomic club.

When you heat any diatomic gas enough, it starts becoming a ‘monatomic’ gas as molecules break down into individual atoms. However, just as a diatomic molecule can break apart into two atoms:

A 2A+A A_2 \to A + A

two atoms can recombine to form a diatomic molecule:

A+AA 2 A + A \to A_2

So in equilibrium, the gas will be a mixture of diatomic and monatomic forms. The exact amount of each will depend on the temperature and pressure, since these affect the likelihood that two colliding atoms stick together, or a diatomic molecule splits apart. The detailed nature of our gas also matters, of course.

But we don’t need to get into these details here! Instead, we can just write down the ‘rate equation’ for the reactions we’re talking about. All the details we’re ignoring will be hiding in some constants called ‘rate constants’. We won’t try to compute these; we’ll leave that to our chemist friends.

A reaction network

To write down our rate equation, we start by drawing a ‘reaction network’. For this, we can be a bit abstract and call the diatomic molecule BB instead of A 2A_2. Then it looks like this:

We could write down the same information using a Petri net:

But today let’s focus on the reaction network! Staring at this picture, we can read off various things:

  • Species. The species are the different kinds of atoms, molecules, etc. In our example the set of species is S={A,B}S = \{A, B\}.

  • Complexes. A complex is a finite sum of species, like AA, or A+AA + A, or for a fancier example, using more efficient notation 2A+3B2 A + 3 B. So, we can think of a complex as a vector v Sv \in \mathbb{R}^S. The complexes that actually show up in our reaction network form a set C SC \subseteq \mathbb{R}^S. In our example, C={A+A,B}C = \{A+A, B\}.

  • Reactions. A reaction is an arrow going from one complex to another. In our example we have two reactions: A+ABA + A \to B and BA+AB \to A + A.

Chemists define a reaction network to be a triple (S,C,T)(S, C, T) where SS is a set of species, CC is the set of complexes that appear in the reactions, and TT is the set of reactions vwv \to w where v,wCv, w \in C. (Stochastic Petri net people call reactions transitions, hence the letter TT.)

So, in our example we have:

  • Species S={A,B}S = \{A,B\}.

  • Complexes: C={A+A,B}C= \{A+A, B\}.

  • Reactions: T={A+AB,BA+A}T = \{A+A\to B, B\to A+A\}.

To get the rate equation, we also need one more piece of information: a rate constant r(τ)r(\tau) for each reaction τT\tau \in \T. This is a nonnegative real number that affects how fast the reaction goes. All the details of how our particular diatomic gas behaves at a given temperature and pressure are packed into these constants!

The rate equation

The rate equation says how the expected numbers of the various species—atoms, molecules and the like—changes with time. This equation is deterministic. It’s a good approximation when the numbers are large and any fluctuations in these numbers are negligible by comparison.

Here’s the general form of the rate equation:

ddtx i= τTr(τ)(n i(τ)m i(τ))x m(τ) \frac{d}{d t} x_i = \sum_{\tau\in T} r(\tau) \, (n_i(\tau)-m_i(\tau)) \, x^{m(\tau)}

Let’s take a closer look. The quantity x ix_i is the expected population of the iith species. So, this equation tells us how that changes. But what about the right hand side? As you might expect, it’s a sum over reactions. And:

  • The term for the reaction τ\tau is proportional to the rate constant r(τ)r(\tau).

  • Each reaction τ\tau goes between two complexes, so we can write it as m(τ)n(τ)m(\tau) \to n(\tau). Among chemists the input m(τ)m(\tau) is called the reactant complex, and the output is called the product complex. The difference n i(τ)m i(τ)n_i(\tau)-m_i(\tau) tells us how many items of species ii get created, minus how many get destroyed. So, it’s the net amount of this species that gets produced by the reaction τ\tau. The term for the reaction τ\tau is proportional to this, too.

  • Finally, the law of mass action says that the rate of a reaction is proportional to the product of the concentrations of the species that enter as inputs. More precisely, if we have a reaction τ\tau with input is the complex m(τ)m(\tau), we define x m(τ)=x 1 m 1(τ)x k m k(τ) x^{m(\tau)} = x_1^{m_1(\tau)} \cdots x_k^{m_k(\tau)}. The law of mass action says the term for the reaction τ\tau is proportional to this, too!

Let’s see what this says for the reaction network we’re studying:

Let’s use x 1(t)x_1(t) (resp. x 2(t)x_2(t)) to stand for the population of species AA (resp. BB). Let the rate constant for the reaction BA+AB \to A + A be r 1r_1, and let the rate constant for A+ABA + A \to B be r 2r_2. Then the rate equation is this:

ddtx 1=2r 1x 22r 2x 1 2 \frac{d}{d t} x_1 = 2 r_1 x_2 - 2 r_2 x_1^2
ddtx 2=r 1x 2+r 2x 1 2 \frac{d}{d t} x_2 = -r_1 x_2 + r_2 x_1^2

This is a bit intimidating. However, we can solve it in closed form thanks to something every physicist knows to treasure: a conserved quantity.

We’ve got two species, AA and BB. But remember, BB is just an abbreviation for a molecule made of two AA atoms. So, the total number of AA atoms is conserved by the reactions A+ABA + A \to B, BA+AB \to A + A. This number is the number of AA‘s plus twice the number of BB’s: x 1+2x 2x_1 + 2x_2. So, this should be a conserved quantity: it should not change with time. Indeed, by adding the first equation above to twice the second, we see:

ddtx 1+2x 2=0 \frac{d}{d t} x_1 + 2x_2 = 0

So, any solution will move along a line

x 1+2x 2=c x_1 + 2 x_2 = c

for some constant cc. We can use this fact to rewrite the rate equation just in terms of x 1x_1:

ddtx 1=2r 1(c2x 1)2r 2x 1 2 \frac{d}{d t} x_1 = 2 r_1 (c - 2x_1) - 2 r_2 x_1^2

and this is a separable differential equation, so we can solve it.

This sort of trick won’t work for more complicated examples. But the idea remains important: the numbers of atoms of various kinds—hydrogen, helium, lithium, and so on—are conserved by chemical reactions, so a solution of the rate equation can’t roam freely in S\mathbb{R}^S. It will be trapped on some hypersurface, which is called the ‘stoichiometric subspace’. And this is very important.

We don’t feel like doing the integral required to solve our rate equation in closed form, because this idea doesn’t generalize too much. On the other hand, we can always solve the rate equation numerically. So let’s try that!

For example, suppose we set r 1=r 2=1r_1 = r_2 = 1. We can plot the solutions for three different choices of initial conditions, say (x 1,x 2)=(0,3),(4,0),(x_1,x_2) = (0,3), (4,0), and (3,3)(3,3). We get these graphs:

It looks like the solution always approaches an equilibrium. We seem to be getting different equilibria for different initial conditions, and the pattern is a bit mysterious. However, something nice happens when we plot the ratio x 1 2/x 2x_1^2 / x_2:

Apparently it always converges to 1. Why should that be? It’s not terribly surprising. With both rate constants equal to 1, the reaction A+ABA + A \to B proceeds at a rate equal to the square of the number of AA‘s, namely x 1 2x_1^2. The reverse reaction proceeds at a rate equal to the number of BB’s, namely x 2x_2. So in equilibrium, we should have x 1 2=x 2x_1^2 = x_2.

But why is the equilibrium stable? In this example we could see that using the closed-form solution, or maybe just common sense. But it also follows from a powerful theorem that handles a lot of reaction networks.

The deficiency zero theorem

It’s called Feinberg’s deficiency zero theorem, and we saw it last time. Very roughly, it says that if our reaction network is ‘weakly reversible’ and has ‘deficiency zero’, the rate equation will have equilibrium solutions that behave about as nicely as you could want.

Let’s see how this works. We need to remember some jargon:

  • Weakly reversible. A reaction network is weakly reversible if for every reaction vwv \to w in the network, there exists a path of reactions in the network starting at ww and leading back to vv.

  • Reversible. A reaction network is reversible if for every reaction vwv \to w in the network, wvw \to v is also a reaction in the network. Any reversible reaction network is weakly reversible. Our example is reversible, since it consists of reactions A+ABA + A \to B, BA+AB \to A + A.

But what about ‘deficiency zero’? We defined that concept last time, but let’s review:

  • Connected component. A reaction network gives a kind of graph with complexes as vertices and reactions as edges. Two complexes lie in the same connected component if we can get from one to the other by a path of reactions, where at each step we’re allowed to go either forward or backward along a reaction. Chemists call a connected component a linkage class. In our example there’s just one:
  • Stoichiometric subspace. The stoichiometric subspace is the subspace Stoch S\mathrm{Stoch} \subseteq \mathbb{R}^S spanned by the vectors of the form wvw - v for all reactions vwv \to w in our reaction network. In our example we have reactions A+ABA + A \to B and BA+AB \to A + A, which give vectors B2AB - 2A and 2AB2A - B, or if you prefer, (2,1)(-2,1) and (2,1)(2,-1). These vectors are linearly dependent, so the stoichiometric subspace has dimension 1.

  • Deficiency. The deficiency of a reaction network is the number of complexes, minus the number of connected components, minus the dimension of the stoichiometric subspace. In our example there are 2 complexes, 1 connected component, and the dimension of the stoichiometric subspace is 1. So, our reaction network has deficiency 2 - 1 - 1 = 0.

So, the deficiency zero theorem applies! What does it say? To understand it, we need a bit more jargon. First of all, a vector x Sx \in \mathbb{R}^S tells us how much we’ve got of each species: the amount of species iSi \in S is the number x ix_i. And then:

  • Stoichiometric compatibility class. Given a vector v Sv\in \mathbb{R}^S, its stoichiometric compatibility class is the subset of all vectors that we could reach using the reactions in our reaction network:
{v+w:wStoch} \{ v + w \; : \; w \in \mathrm{Stoch} \}

In our example, where the stoichiometric subspace is spanned by (2,1)(2,-1), the stoichiometric compatibility class of the vector (a,b)(a,b) is the line consisting of points

(x 1,x 2)=(a,b)+s(2,1) (x_1, x_2) = (a,b) + s(2,-1)

where the parameter ss ranges over all real numbers. Notice that this line can also be written as

x 1+2x 2=c x_1 + 2x_2 = c

We’ve already seen that if we start with initial conditions on such a line, the solution will stay on this line. And that’s how it always works: as time passes, any solution of the rate equation stays in the same stoichiometric compatibility class!

In other words: the stoichiometric subspace is defined by a bunch of linear equations, one for each linear conservation law that all the reactions in our network obey. Here by a linear conservation law we mean a law saying that some linear combination of the numbers of species does not change.

Next:

  • Positivity. A vector in S\mathbb{R}^S is positive if all its components are positive; this describes a a container of chemicals where all the species are actually present. The positive stoichiometric compatibility class of x Sx\in \mathbb{R}^S consists of all positive vectors in its stoichiometric compatibility class.

We now finally have enough jargon in our arsenal to state this result:

Zero Deficiency Theorem (Feinberg). If a reaction network is weakly reversible and the rate constants are positive, the rate equation has exactly one equilibrium solution in each positive stoichiometric compatibility class. Any sufficiently nearby solution that starts in the same class will approach this equilibrium as t+t \to +\infty.

In our example, this theorem says there’s just one positive equilibrium (x 1,x 2)(x_1,x_2) in each line

x 1+2x 2=c x_1 + 2x_2 = c

We can find it by setting the time derivatives to zero:

ddtx 1=2r 1x 22r 2x 1 2=0 \frac{d}{d t}x_1 = 2 r_1 x_2 - 2 r_2 x_1^2 = 0
ddtx 2=r 1x 2+r 2x 1 2=0 \frac{d}{d t} x_2 = -r_1 x_2 + r_2 x_1^2 = 0

Solving these, we get

x 1 2x 2=r 1r 2 \frac{x_1^2}{x_2} = \frac{r_1}{r_2}

So, these are our equilibrium solutions. It’s easy to verify that indeed, there’s one of these in each stoichiometric compatibility class x 1+2x 2=cx_1 + 2x_2 = c. And the zero deficiency theorem also tells us that any sufficiently nearby solution that starts in the same class will approach this equilibrium as tt \to \infty.

This partially explains what we saw before in our graphs. It shows that in the case r 1=r 2=1r_1 = r_2 = 1, any solution that starts by nearly having

x 1 2x 2=1 \frac{x_1^2}{x_2} = 1

will actually have

lim t+x 1 2x 2=1 \lim_{t \to +\infty} \frac{x_1^2}{x_2} = 1

But in fact, in this example we don’t even need to start near the equilibrium for our solution to approach the equilibrium! What about in general? We don’t know, but just to get the ball rolling, we’ll risk the following wild guess:

Conjecture. If a reaction network is weakly reversible and the rate constants are positive, the rate equation has exactly one equilibrium solution in each positive stoichiometric compatibility class, and any positive solution that starts in the same class will approach this equilibrium as t+t \to +\infty.

If anyone knows a proof or counterexample, we’d be interested. If this result were true, it would really clarify the dynamics of reaction networks in the zero deficiency case.

category: blog