Network Theory (Part 18)

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, the deficiency zero theorem, which 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_2 \to A + A$

two atoms can recombine to form a diatomic molecule:

$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 $B$ instead of $A_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\}$ .

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

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

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

So, in our example we have:

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

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

• Reactions: $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(\tau)$ for each reaction $\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:

$\displaystyle{ \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_i$ is the expected population of the $i$ th 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(\tau).$

• Each reaction $\tau$ goes between two complexes, so we can write it as $m(\tau) \to n(\tau).$ Among chemists the input $m(\tau)$ is called the reactant complex, and the output is called the product complex. The difference $n_i(\tau)-m_i(\tau)$ tells us how many items of species $i$ 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(\tau)$ , we define $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 write $x_1(t)$ for the number of $A$ atoms and $x_2(t)$ for the number of $B$ molecules. Let the rate constant for the reaction $B \to A + A$ be $\alpha,$ and let the rate constant for $A + A \to B$ be $\beta,$ Then the rate equation is this:

$\displaystyle{\frac{d}{d t} x_1 = 2 \alpha x_2 - 2 \beta x_1^2 }$

$\displaystyle{\frac{d}{d t} x_2 = -\alpha x_2 + \beta x_1^2 }$

This is a bit intimidating. However, we can solve it in closed form thanks to something very precious: a conserved quantity.

We’ve got two species, $A$ and $B$ . But remember, $B$ is just an abbreviation for a molecule made of two $A$ atoms. So, the total number of $A$ atoms is conserved by the reactions in our network. This is the number of A‘s plus twice the number of B‘s: $x_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:

$\displaystyle{\frac{d}{d t} \left( x_1 + 2x_2 \right) = 0 }$

As a consequence, any solution will stay on a line

$x_1 + 2 x_2 = c$

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

$\displaystyle{ \frac{d}{d t} x_1 = \alpha (2c - x_1) - 2 \beta x_1^2 }$

This is a separable differential equation, so we can solve it if we can figure out how to do this integral

$\displaystyle{ t = \int \frac{d x_1}{\alpha (2c - x_1) - 2 \beta x_1^2 } }$

and then solve for $x_1.$

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 $\mathbb{R}^S.$ It will be trapped on some hypersurface, which is called a ‘stoichiometric compatibility class’. 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 $\alpha = \beta = 1.$ We can plot the solutions for three different choices of initial conditions, say $(x_1,x_2) = (0,3), (4,0),$ and $(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_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 + A \to B$ proceeds at a rate equal to the square of the number of $A$ ‘s, namely $x_1^2$ . The reverse reaction proceeds at a rate equal to the number of $B$ ‘s, namely $x_2.$ So in equilibrium, we should have $x_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 $v \to w$ in the network, there exists a path of reactions in the network starting at $w$ and leading back to $v.$

• Reversible. A reaction network is reversible if for every reaction $v \to w$ in the network, $w \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 + A \to B,$ $B \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 $\mathrm{Stoch} \subseteq \mathbb{R}^S$ spanned by the vectors of the form $w - v$ for all reactions $v \to w$ in our reaction network. This subspace describes the directions in which a solution of the rate equation can move. In our example, it’s spanned by $B - 2 A$ and $2 A - B$ , or if you prefer, $(-2,1)$ and $(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 \in \mathbb{R}^S$ tells us how much we’ve got of each species: the amount of species $i \in S$ is the number $x_i.$ And then:

• Stoichiometric compatibility class. Given a vector $v\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 \; : \; w \in \mathrm{Stoch} \}$

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

$(x_1, x_2) = (v_1,v_2) + s(2,-1)$

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

$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 a linear conservation law is a law saying that some linear combination of the numbers of species does not change.

• Positivity. A vector in $\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\in \mathbb{R}^S$ consists of all positive vectors in its stoichiometric compatibility class.

We finally have enough jargon in our arsenal to state the zero deficiency theorem. We’ll only state the part we need today:

Zero Deficiency Theorem (Feinberg). If a reaction network has deficiency zero and 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 \to +\infty.$

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

$x_1 + 2x_2 = c$

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

$\displaystyle{ \frac{d}{d t}x_1 = 2 \alpha x_2 - 2 \beta x_1^2 = 0 }$

$\displaystyle{ \frac{d}{d t} x_2 = -\alpha x_2 + \beta x_1^2 = 0 }$

Solving these, we get

$\displaystyle{ \frac{x_1^2}{x_2} = \frac{\alpha}{\beta} }$

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 = 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 $t \to \infty$ .

This partially explains what we saw before in our graphs. It shows that in the case $\alpha = \beta = 1,$ any solution that starts by nearly having

$\displaystyle{\frac{x_1^2}{x_2} = 1}$

will actually have

$\displaystyle{\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 has deficiency zero and 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 \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.

Next time we’ll talk about this same reaction network from a stochastic point of view, where we think of the atoms and molecules as reacting in a probabilistic way. And we’ll see how the conservation laws we’ve been talking about today are related to Noether’s theorem for Markov processes!

This entry was posted on Monday, November 21st, 2011 at 7:07 am and is filed under chemistry, mathematics, networks, physics. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.

27 Responses to Network Theory (Part 18)

David Corfield says:

21 November, 2011 at 12:20 pm

In your conjecture, presumably you’re still taking deficiency to be zero. It would be nice to see the simplest deficiency one network where there fails to be a single fixed point per coset. Feinberg wrote about deficiency one networks, as mentioned by Gunawardena. That seems a very good paper, if one had time. That idea that his $\Psi$ is “trying to be linear” looks important.

I guess what your conjecture needs is prevention of states heading off to the boundary of the non-negative cone, where one or more of the concentrations tends to zero.

Reply
- John Baez says:
  
  21 November, 2011 at 12:49 pm
  
  Yikes, we left the ‘deficiency zero’ hypothesis out of our statement of the deficiency zero theorem, and also our conjecture! Thanks for catching that—I’ll put it in.
  
  I’ve been sadly slow understanding the proof of the deficiency zero theorem or Feinberg’s work on deficiency one networks. Gunawardena’s paper is very good, but Feinberg’s original papers also look very good (and free):
  
  • Martin Feinberg, Chemical reaction network structure and the stability of complex isothermal reactors: I. The deficiency zero and deficiency one theorems, Chemical Engineering Science 42 (1987), 2229–2268.
  
  • Martin Feinberg, Lectures on reaction networks, 1979.
  
  I’ve just been slow about reading them!
  
  I guess what your conjecture needs is prevention of states heading off to the boundary of the non-negative cone, where one or more of the concentrations tends to zero.
  
  That’s part of it: that’s what can easily go wrong if we drop the ‘weakly reversible’ condition. Periodic orbits are already ruled out. But there’s also the possibility, in principle, of more complicated chaotic motions!
  
  A standard trick in situations like this would be to cook up a ‘Lyapunov function’ on each stoichiometric compatibility class, in order to show each class contains a globally asymptotically stable equilibrium—which just means an equilibrium that all other solutions (in that class) approach as $t \to +\infty$ . I should look and see what Feinberg did, and whether or why he didn’t or couldn’t do that.
  
  Reply
- John Baez says:
  
  21 November, 2011 at 1:17 pm
  
  David wrote:
  
  It would be nice to see the simplest deficiency one network where there fails to be a single fixed point per coset.
  
  Let’s see. Feinberg mentions a nice example on page 5 here: it’s called the Edelstein network:
  
  $A \leftrightarrow A + A$
  $A + B \leftrightarrow C \leftrightarrow B$
  
  This has 5 complexes and 2 components, and the stoichiometric subspace has dimension 2, since it’s the span of the vectors $A$ and $B - C$ . So, the deficiency of this reaction network is 5 – 2 – 2 = 1.
  
  He says that for certain choices of the rate constants, the set of equilibrium is a roughly S-shaped curve that hits one of the stoichiometric compatibility classes—which are planes—in three points! For example, he says these rate constants work:
  
  $8.5: A \to A + A$
  $1: A + A \to A$
  $1: A + B \to C$
  $1: C \to A + B$
  $1 : C \to B$
  $0.2: B \to C$
  
  I’d like to play with this and other examples myself, soon. He lists some more. These ‘play’ networks, as he calls them, all seem unrealistic in one way or another—at least in chemistry. In particular, $A \to A + A$ may be fine for amoebas, I don’t think it happens in chemistry—unless, just as our amoebas are supplied with food, our chemical is a self-replicating molecule floating in a ‘bath’ of raw materials that we don’t take explicit account of here.
  
  Reply
- Blake Stacey says:
  
  21 November, 2011 at 1:28 pm
  
  Periodic orbits are already ruled out.
  
  I clearly don’t yet understand this stuff well enough, because the ruling-out of periodic orbits wasn’t transparently obvious. (The first kind of counterexample I’d try to find to the conjecture would be a reaction network whose dynamics go to a limit cycle of some sort.)
  
  But there’s also the possibility, in principle, of more complicated chaotic motions!
  
  Time to dust off that part of my brain which learned about the Belousov–Zhabotinsky reaction….
  
  Reply
  - John Baez says:
    
    21 November, 2011 at 1:35 pm
    
    Blake wrote:
    
    I clearly don’t yet understand this stuff well enough, because the ruling-out of periodic orbits wasn’t transparently obvious.
    
    Just so nobody gets confused: periodic orbits are ruled out for deficiency zero, weakly reversible reaction networks—this is part of the deficiency zero theorem as stated last time.
    
    This theorem is not obvious to me, so we’re in the same boat: I need to understand what’s so great about ‘deficiency zero’.
    
    Reaction networks violating these hypotheses can have periodic orbits:
    
    • Daniel B. Forger, Signal processing in cellular clocks, PNAS 108 (2011), 4281-4285.
    
    Time to dust off that part of my brain which learned about the Belousov–Zhabotinsky reaction…
    
    Me too! It’s not chaotic, but you might like the reaction on page 10 here, which is based on the ‘Brusselator’. It has an unstable equilibrium inside a stable limit cycle.
    
    (I used to see ‘Brusselator’ and ‘Belousov–Zhabotinsky’ in the same books and papers, maybe just because they’re both examples of interestingly complex chemical reactions. But somehow ‘Brusselator’ sounds like a science-fiction weapon made up by a kid who hates Brussels sprouts.)
    
    Reply
  - David Corfield says:
    
    21 November, 2011 at 2:50 pm
    
    I take it those chemical oscillators aren’t weakly reversible. It doesn’t look like it from here.
    
    Actually, that’s really a very strong condition, if you think about it, that every reaction between complexes may be reversed along some path. I see Gunawardena points out that
    
    …there are more intricate forms of reversibility in which both complexes and species play a part. (p. 16)
    
    He refers to P. M. Schlosser’s thesis, A graphical determination of the possibility of multiple steady states in complex isothermal CFSTRs. PhD thesis, University of Rochester, 1988.
    
    Reply
  - Blake Stacey says:
    
    21 November, 2011 at 10:50 pm
    
    To quote the song, I don’t know much about the BZ reaction, but I’ve heard it can display chaotic behavior when it’s pushed into the right regime in a machine called a “continuous-flow, stirred-tank reactor”. Some historical references pulled from Strogatz’s Nonlinear Dynamics and Chaos (1994) are:
    
    R. A. Schmitz, K. R. Graziani, and J. L. Hudson (1977), “Experimental evidence of chaotic states in the Belousov–Zhabotinskii reaction” J. Chem. Phys. 67, 3040, doi:10.1063/1.435267.
    
    R. H. Simoyi, A. Wolf and H. L. Swinney (1982), “One-dimensional dynamics in a multicomponent chemical reaction” Physical Review Letters 49, 245, doi:10.1103/PhysRevLett.49.245.
    
    J. C. Roux, R. H. Simoyi and H. L. Swinney (1983), “Observation of a strange attractor” Physica D 8, 257, doi:10.1016/0167-2789(83)90323-8.
    
    The “Brusselator” is a simplified model which tries to capture the BZ dynamics. The Oregonator, invented more recently, appears to be a more chemically realistic version.
    
    Reply
  - Blake Stacey says:
    
    21 November, 2011 at 11:00 pm
    
    The rate equations for the Oregonator are as follows:
    
    $dX/dt = k_1 A Y - k_2 X Y + k_3 A X - 2 k_4 X^2$
    
    $dY/dt = -k_1 A Y - k_2 X Y + 1/2 k_5 f B Z$
    
    $dZ/dt = 2k_3 A X - k_5 B Z$
    
    Here, $X$ denotes the concentration of HBrO2, $Y$ that of the negative bromide ion and $Z$ that of the positive cerium ion. The other letters stand for concentrations which change more slowly and are often taken to be constant.
    
    (Not a chemist — I’m just reading all this out of Scholarpedia.)
    
    Reply
  - John Baez says:
    
    22 November, 2011 at 12:04 am
    
    Here is the stable limit cycle of the chemical reaction network Blake pointed out, the ‘Oregonator’:
    
    This is taken from
    
    • Richard J. Field, Oregonator, Scholarpedia.
    
    Reply
  - John Baez says:
    
    22 November, 2011 at 12:35 am
    
    David wrote:
    
    I take it those chemical oscillators aren’t weakly reversible. It doesn’t look like it from here.
    
    Thanks for that link! Yes, the way they describe the Bray–Liebhafsky reaction, there are reactions that turn hydrogen peroxide to water but not vice versa:
    
    5 H₂O2 + I₂ → 2 IO^3- + 2 H⁺ + 4 H₂O
    
    5 H₂O₂ + 2 IO^3- + 2 H⁺ → I₂ + 5 O₂ + 6 H₂O
    
    Here the decay of hydrogen peroxide to water is powering a cycle where iodine gets oxidized to iodate and then iodate gets reduced back to iodine.
    
    So, it’s a bit similar to how the wheels on a cart turn around as the cart rolls downhill. They only turn as it’s rolling; when the system reaches equilibrium at the bottom of the hill, they stop.
    
    So, I think the Bray–Liebhafsky reaction will only oscillate as long as we have plenty of hydrogen peroxide around; in the $t \to +\infty$ limit it all becomes water, and things become boring.
    
    But as John Maynard Keynes noted, “in the long run we’re all dead”. So, we shouldn’t get too focused on $t \to +\infty$ behavior if we’re studying biology.
    
    Actually, that’s really a very strong condition, if you think about it, that every reaction between complexes may be reversed along some path.
    
    You’re right: so-called ‘weak reversibility’ is pretty strong.
    
    For example, this reaction network is not weakly reversible:
    
    A → B
    
    A + B → A + A
    
    However, as long as we have at least one A around to serve as a ‘catalyst’, we can use these reactions to turn any number of A’s into the same number of B’s, or vice versa.
    
    On the other hand: in real life all chemical reaction networks are not just weakly reversible but actually reversible. If a chemical reaction is possible, so is its reverse! But if the rate constant for the reverse reaction is sufficiently low, we often neglect it.
    
    Reply
    - Anonymous says:
      
      2 December, 2011 at 7:26 am
      
      John,
      
      Given your interest in leaves over hypercubes (e.g., 1st network th post), perhaps the Belousov–Zhabotinsky reaction would stimulate some ideas? It is an example of a class of reactions that under certain circumstances, exhibit fundamental self organizing properties (like leaves!).
    - John Baez says:
      
      2 December, 2011 at 8:07 am
      
      Yes, I would like to think about this! I’ve started with a lot of general formalism concerning chemical reaction networks, and some examples that are fairly dull in the sense that the solutions converge to an equilibrium, but I want to start studying self-organizing systems… as well as expand the ‘network theory’ ideas to cover other subjects. I’ve done more on the latter, so far. There’s a lot of stuff already done that I need to explain. But thanks for pushing me towards the former!
Uwe Stroinski says:

21 November, 2011 at 6:51 pm

Concerning your conjecture the following might be helpful:

• David F. Anderson, Global asymptotic stability for a class of nonlinear chemical equations, SIAM J. Appl. Math. 68(5) (May 2008), 1464 – 1476.

He generalizes a result of E. Sontag (from 2001) that ensures global asymptotic stability if there are no equilibria on the boundary of the positive compatibility classes.

Reply
- John Baez says:
  
  22 November, 2011 at 12:47 pm
  
  Thanks a million, Uwe! I added a link to the arXiv version.
  
  In response to a question by David, I guessed that the key to proving stability would be coming up with a suitable Lyapunov function. That’s a pretty obvious guess, of course: no stroke of genius required there. But I’m happy to see that in the paper you cite, a specific candidate Lyapunov function is written down in equation (3)—and David will be just as pleased as me to see that the formula is reminiscent of entropy!
  
  Maybe it’s more like ‘free energy’, but no matter: we’ve already seen that certain simple reaction networks show up in evolutionary game theory, and there the Lyapunov function is indeed a close relative of entropy: namely, the information the species has left to learn about its environment! This keeps dropping with the passage of time, and the species is in equilibrium when it has nothing left to learn.
  
  As we saw back then, Marc Harper has a nice paper on this subject. The most relevant section is the one called ““Kullback-Leibler Divergence is a Lyapunov function for the Replicator Dynamic”.
  
  I wrote about this in a post on information geometry. And I’m happy now because my interest in information geometry and my interest in reaction networks are starting to overlap!
  
  Reply
Manoj Gopalkrishnan says:

22 November, 2011 at 4:01 am

Dear John,

I was going to point you to the global attractor conjecture, but I see Uwe Stroinski has already done this. Your conjecture, and some generalizations of it, are well-known open questions in the community. In fact, in their 1972 paper, Horn and Jackson claimed a proof for a similar conjecture, but Horn found a problem with their proof, and withdrew it in 1974, and the problem has remained open since then. A generalization of the conjecture was posed in 1989 by Marty Feinberg. This has come to be called the “persistence conjecture.”

The problem is, in essence, the following: starting from a strictly positive vector of concentrations, can species become extinct? If they can not, one says the network is persistent. Feinberg’s conjecture is that if the network is weakly-reversible, then it must be persistent: species can not go extinct. It is not hard to show that this implies the global attractor conjecture, not just for deficiency zero, but for all complex-balanced systems (roughly, systems that admit a “free energy” function).

Here’s a far-from-comprehensive summary of the state of the art: even the three-dimensional version (the stoichiometric compatibility class is 3-dimensional) of the persistence conjecture is open. The two-dimensional version has been settled due to recent work by Craciun, Nazarov, and Pantea http://arxiv.org/abs/1010.3050 . David Anderson also has a very nice recent paper http://arxiv.org/abs/1101.0761 that settles a special case (the graph is strongly-connected) of this result. I can prove (http://arxiv.org/abs/1006.3627 ) the persistence conjecture in the absence of catalysis in the network.

Reply
- John Baez says:
  
  22 November, 2011 at 1:24 pm
  
  Thanks for all the info, Manoj! This is really great. I hadn’t expected my conjecture was new, so I’m happy to see it might be true, and intrigued that it’s a well-known open puzzle.
  
  This is very interesting:
  
  Feinberg’s conjecture is that if the network is weakly-reversible, then it must be persistent: species can not go extinct. It is not hard to show that this implies the global attractor conjecture, not just for deficiency zero, but for all complex-balanced systems (roughly, systems that admit a “free energy” function).
  
  This is basically what David guessed:
  
  I guess what your conjecture needs is prevention of states heading off to the boundary of the non-negative cone, where one or more of the concentrations tends to zero.
  
  So, 10 points for David! But I’ll need to read why it’s actually true, because a priori I could imagine problems other than one or more concentrations going to zero (= species going extinct).
  
  Anyway, I’ve got my reading cut out for me!
  
  Reply
  - Manoj Gopalkrishnan says:
    
    22 November, 2011 at 3:11 pm
    
    It’s not that hard to see. Here’s the basic idea. Consider the omega-limit set of a trajectory. Note the following two facts about this set:
    
    1. Every point on this omega-limit set has the same value of the Lyapunov function.
    
    2. The omega-limit set is invariant for the dynamics.
    
    Therefore, the time derivative of the Lyapunov function is zero for points in the omega-limit set. Since we have a strict Lyapunov function in the relative interior of the invariant polyhedron, this omega-limit set can not contain any positive point, except the isolated positive equilibrium points. So if there are no omega-limit points on the boundary, we’re done.
    
    Reply
    - Manoj Gopalkrishnan says:
      
      22 November, 2011 at 3:13 pm
      
      I should add that trajectories can not escape to infinity because the level sets of our Lyapunov function are bounded.
Manoj Gopalkrishnan says:

22 November, 2011 at 8:45 am

By the way, I really hate saying “positive stoichiometric compatibility class.” I prefer the terms “conservation class” instead of stoichiometric compatibility class (the affine space), and “invariant polyhedron” for the intersection of the conservation class with the positive orthant (possibly very close to what you call the positive stoichiometric compatibility class, I’ll have to check your definition).

Reply
- John Baez says:
  
  22 November, 2011 at 1:11 pm
  
  Yes, I have been using ‘positive stoichiometric compatibility class’ to mean the intersection of a stochiometric compatibility class (an affine space $x + \mathrm{Stoch}$ , in my notation) with the open positive orthant (the set $(0,\infty)^S$ , in my notation). I agree that this is long and clunky, but I’ve seen various experts do it and I haven’t dared deviate yet.
  
  I agree that mathematicians would find ‘invariant polyhedron’ vastly more appealing. For one thing, few of them know what ‘stoichiometric’ means, and it’s a bit of a jawbreaker. For another, ‘invariant polyhedron’ will instantly make some sense to them.
  
  To go along with ‘invariant polyhedron’, one might say ‘invariant coset’ for a stochiometric compatibility class, since it’s a coset of the subspace $\mathrm{Stoch}$ in the vector space $\mathbb{R}^S$ . Then it all sounds very much like math.
  
  Reply
Weekly Picks « Mathblogging.org — the Blog says:

1 December, 2011 at 4:10 am

[…] in the ongoing series about reaction networks, Azimuth had a rather self-contained post about the practical applications of the deficiency zero theorem […]

Reply
Network Theory (Part 19) « Azimuth says:

18 July, 2012 at 9:56 am

[…] Last time we started looking at a simple example: a diatomic gas […]

Reply
Arjun Jain says:

4 June, 2013 at 9:18 am

If s is the number of species, and r the dimensionality of $\mathrm{Stoch};$ in the language of linear algebra, $\mathrm{Stoch}$ can be identified with the row space. Then as the null space is the orthogonal complement of the row space in $\mathbb{R}^s,$ we have s-r independent null space vectors, and therefore s-r linear conservation laws. So, in deficiency 0 networks the number of connected components is equal to the number of conservation laws. Is this correct?

Reply
- John Baez says:
  
  4 June, 2013 at 4:06 pm
  
  The number of linear conservation laws is the number of species minus the dimension of the stoichiometric subspace. Even better, the procedure given here for determining the stoichiometric subspace can be used to find all the linear conservation laws. So, this answers your earlier question about how to find conservation laws.
  
  So, in deficiency 0 networks the number of connected components is equal to the number of conservation laws. Is this correct?
  
  Let’s see. The deficiency is the number of complexes minus the number of connected components minus the dimension of the stochiometric subspace. So, if it’s zero,
  
  • The number of complexes equals the number of connected components plus the dimension of the stochiometric subspace.
  
  On the other hand,
  
  • The number of linear conservation laws is the number of species minus the dimension of the stoichiometric subspace.
  
  Putting these equations together, I can get an equation involving the number of connected components and the number of conservation laws. But it doesn’t seem to be the equation you suggest! What equation do we actually get?
  
  Reply
Arjun Jain says:

4 June, 2013 at 4:59 pm

Yes, I got a little mixed up.
We get : no. of complexes – no. of species = no. of connected components – no. of conservation laws.

Reply
- John Baez says:
  
  4 June, 2013 at 7:42 pm
  
  Okay, that looks great! Think about some examples. And think about some counterexamples when the reaction network doesn’t have deficiency zero. I think it should make intuitive sense after a while.
  
  Reply
Quantum Theory and Complex Networks – QuamPlexity.org says:

24 April, 2016 at 6:03 pm

[…] An example of the deficiency zero theorem: a diatomic gas, Network Theory Part 18 […]

Reply

You can use Markdown or HTML in your comments. You can also use LaTeX, like this: $latex E = m c^2 $. The word 'latex' comes right after the first dollar sign, with a space after it. Cancel reply

Enter your comment here...

Fill in your details below or click an icon to log in:

Email (required) (Address never made public)

Name (required)

Website

You are commenting using your WordPress.com account. ( Log Out / Change )

You are commenting using your Twitter account. ( Log Out / Change )

You are commenting using your Facebook account. ( Log Out / Change )

Cancel

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.

	Jim Stuttard on Azimuth Project News
	Toby Bartels on Chemistry and Invariant T…
	Wolfgang on Chemistry and Invariant T…
	Raphael on Chemistry and Invariant T…
	Bill Bottenberg on Chemistry and Invariant T…
	Wyrd Smythe on Four ‘Universes’
	Phillip Helbig on Four ‘Universes’
	Phillip Helbig on Four ‘Universes’
	Phillip Helbig on Four ‘Universes’
	Wyrd Smythe on The Galactic Center

Azimuth