The summer before last, I invited Brendan Fong to Singapore to work with me on my new ‘network theory’ project. He quickly came up with a nice new proof of a result about mathematical chemistry. We blogged about it, and I added it to my book, but then he became a grad student at Oxford and got distracted by other kinds of networks—namely, Bayesian networks.
So, we’ve just now finally written up this result as a self-contained paper:
• John Baez and Brendan Fong, Quantum techniques for studying equilibria in chemical reaction networks.
Check it out and let us know if you spot mistakes or stuff that’s not clear!
The idea, in brief, is to use math from quantum field theory to give a somewhat new proof of the Anderson–Craciun–Kurtz theorem.
This remarkable result says that in many cases, we can start with an equilibrium solution of the ‘rate equation’ which describes the behavior of chemical reactions in a deterministic way in the limit of a large numbers of molecules, and get an equilibrium solution of the ‘master equation’ which describes chemical reactions probabilistically for any number of molecules.
The trick, in our approach, is to start with a chemical reaction network, which is something like this:
and use it to write down a Hamiltonian describing the time evolution of the probability that you have various numbers of each kind of molecule: A, B, C, D, E, … Using ideas from quantum mechanics, we can write this Hamiltonian in terms of annihilation and creation operators—even though our problem involves probability theory, not quantum mechanics! Then we can write down the equilibrium solution as a ‘coherent state’. In quantum mechanics, that’s a quantum state that approximates a classical one as well as possible.
All this is part of a larger plan to take tricks from quantum mechanics and apply them to ‘stochastic mechanics’, simply by working with real numbers representing probabilities instead of complex numbers representing amplitudes!
I should add that Brendan’s work on Bayesian networks is also very cool, and I plan to talk about it here and even work it into the grand network theory project I have in mind. But this may take quite a long time, so for now you should read his paper:
• Brendan Fong, Causal theories: a categorical perspective on Bayesian networks.
Posted by John Baez 
is a finite set and 
is a probability distribution on that set. Suppose
Then the expected value of 
with respect to 
of the function 
is bigger. We pronounce 
like this: “the expected value of 
If you get heads, I’ll give you 
more points. If you get tails, I’ll take away 
, since I’m multiplying your score by 6. If you lose, your payoff will be 
So, the expected value of your payoff is still zero:
But it’s much bigger than the number of galaxies in the observable universe, which is about 
, or even the number of stars in the observable universe, which is roughly 
It’s impressive that we can hold such a big number in our hand… in the form of possible ways to order a deck of cards!
of them, the probability that they’re in any particular order is 
cards from a well-shuffled deck of 
cards. You’ll get a subset of cards—though card players call this a hand of cards instead of a subset. And, there are 
Just multiply the top and bottom of that big fraction by 
of landing heads up and 
of landing tails up;
of landing heads up and 
of landing tails up;
of landing heads up and 
of landing tails up. 
is
? Because it’ll help us see the pattern for doing harder problems like this!
possible events that can happen. Each of these is equally likely if it’s a fair coin and the flips are independent. So each has probability 
by this. So the answer is
is the number of ways of choosing a subset of 
Remember that 
, the Cartesian product of the sets 
, is the set of all ordered pairs 
where 
and 
:
and 
is a probability distribution on a set of events 
on 
instead of 
on 
with
on 
This tells us the answer to our question. We can work it out:
while the second is described by 
How likely is it that the first lands heads up and the second lands tails up? We can answer questions like this if the coin flips are independent. We do this by multiplying 
on 
can be either 
or 
we’ll call its probability 
called a probability to each event 
such that:
we have
Then by rule 1) above,
will occur. We define this to be
and 
must be equal: 
![\mathrm{rate}(T) = u [A] [B]](http://s0.wp.com/latex.php?latex=%5Cmathrm%7Brate%7D%28T%29+%3D+u+%5BA%5D+%5BB%5D+&bg=ffffff&fg=333333&s=0 "\mathrm{rate}(T) = u [A] [B]")
![\mathrm{rate}(T') = v [C] [D]](http://s0.wp.com/latex.php?latex=%5Cmathrm%7Brate%7D%28T%27%29+%3D+v+%5BC%5D+%5BD%5D+&bg=ffffff&fg=333333&s=0 "\mathrm{rate}(T') = v [C] [D]")
![[X]](http://s0.wp.com/latex.php?latex=%5BX%5D&bg=ffffff&fg=333333&s=0 "[X]")
means the concentration of 
![u [A] [B] = v [C] [D]](http://s0.wp.com/latex.php?latex=u+%5BA%5D+%5BB%5D+%3D+v+%5BC%5D+%5BD%5D+&bg=ffffff&fg=333333&s=0 "u [A] [B] = v [C] [D]")
![\displaystyle{ \frac{[A][B]}{[C][D]} = \frac{v}{u} = K }](http://s0.wp.com/latex.php?latex=%5Cdisplaystyle%7B+%5Cfrac%7B%5BA%5D%5BB%5D%7D%7B%5BC%5D%5BD%5D%7D+%3D+%5Cfrac%7Bv%7D%7Bu%7D+%3D+K+%7D&bg=ffffff&fg=333333&s=0 "\displaystyle{ \frac{[A][B]}{[C][D]} = \frac{v}{u} = K }")
is the equilibrium constant for the reaction pair. 
![u [A]^2 = v [D]](http://s0.wp.com/latex.php?latex=u+%5BA%5D%5E2+%3D+v+%5BD%5D+&bg=ffffff&fg=333333&s=0 "u [A]^2 = v [D]")
and 
be the population counts for A and B, and let
is an invariant property of the system. The reactions cannot change it, because they are just shuffling the units of A from one arrangement to the other. By way of contrast, 
, we get:![[N] = [A] + 2 [D]](http://s0.wp.com/latex.php?latex=%5BN%5D+%3D+%5BA%5D+%2B+2+%5BD%5D+&bg=ffffff&fg=333333&s=0 "[N] = [A] + 2 [D]")
![[N]](http://s0.wp.com/latex.php?latex=%5BN%5D&bg=ffffff&fg=333333&s=0 "[N]")
is the concentration of the units of A.
and 
, we can then solve for the concentrations at equilibrium:![\displaystyle{u [A]^2 = v [D] = v ([N] - [A]) / 2 }](http://s0.wp.com/latex.php?latex=%5Cdisplaystyle%7Bu+%5BA%5D%5E2+%3D+v+%5BD%5D+%3D+v+%28%5BN%5D+-+%5BA%5D%29+%2F+2+%7D&bg=ffffff&fg=333333&s=0 "\displaystyle{u [A]^2 = v [D] = v ([N] - [A]) / 2 }")
![\displaystyle{2 u [A]^2 + v [A] - v [N] = 0 }](http://s0.wp.com/latex.php?latex=%5Cdisplaystyle%7B2+u+%5BA%5D%5E2+%2B+v+%5BA%5D+-+v+%5BN%5D+%3D+0+%7D&bg=ffffff&fg=333333&s=0 "\displaystyle{2 u [A]^2 + v [A] - v [N] = 0 }")
![\displaystyle{[A] = (-v \pm \sqrt{v^2 + 8 u v [N]}) / 4 u }](http://s0.wp.com/latex.php?latex=%5Cdisplaystyle%7B%5BA%5D+%3D+%28-v+%5Cpm+%5Csqrt%7Bv%5E2+%2B+8+u+v+%5BN%5D%7D%29+%2F+4+u+%7D&bg=ffffff&fg=333333&s=0 "\displaystyle{[A] = (-v \pm \sqrt{v^2 + 8 u v [N]}) / 4 u }")
![[A]](http://s0.wp.com/latex.php?latex=%5BA%5D&bg=ffffff&fg=333333&s=0 "[A]")
can’t be negative, only the positive square root is valid.
:![\displaystyle{[A] = (\sqrt{8 [N] + 1} - 1) / 4 }](http://s0.wp.com/latex.php?latex=%5Cdisplaystyle%7B%5BA%5D+%3D+%28%5Csqrt%7B8+%5BN%5D+%2B+1%7D+-+1%29+%2F+4+%7D&bg=ffffff&fg=333333&s=0 "\displaystyle{[A] = (\sqrt{8 [N] + 1} - 1) / 4 }")
![\displaystyle{[D] = ([N] - [A]) / 2 }](http://s0.wp.com/latex.php?latex=%5Cdisplaystyle%7B%5BD%5D+%3D+%28%5BN%5D+-+%5BA%5D%29+%2F+2+%7D&bg=ffffff&fg=333333&s=0 "\displaystyle{[D] = ([N] - [A]) / 2 }")
like the trajectory of a thrown stone. But there are situations where we have some kind of noise, a chaotic, fluctuating influence, that we would like to take into account. This could be, for example, turbulence in the air flow around a rocket. Or, in our case, short term fluctuations of the weather of the earth. If we take this into account, we get an equation of the form
is a ‘random function’ which models the noise. Typically this noise is just a simplified way to take into account rapidly changing fine-grained aspects of the system at hand. This way we do not have to explicitly model these aspects, which is often impossible. 
is 
is a bit like a function, but often distributions are too nasty to have well-defined values at points! Instead, all that makes sense are expressions like this:
and then integrate it. Indeed, we specify a distribution just by saying what all these integrals equal. For example, the 
is a distribution defined by
, but its integral should equal 1. So, if you try to graph it, the region under the graph should be an infinitely skinny infinitely tall spike whose area is 1:
is Brownian noise, so is 
for all 
.
is some fixed number, our ‘time step’. Then we can define
are independent Gaussian random variables with mean zero and standard deviation 
if we are given some initial value 
. And if the time step 
is small enough, we can hope to get a ‘good approximation to the true solution’. Of course defining what we mean by a ‘good approximation’ is tricky here… but I think it’s more important to just plunge in and see what happens, to get a feel for what’s going on here.
and 
are constants we get to choose. If we set them both to zero we get:
and an unstable equilibrium in between at 
. So, this is a bistable model similar to the one we’ve been studying, but mathematically simpler! 
, while the red curve is 
. Here is a simulation with a low level of noise:
to 
, or vice versa. The solution 
and the noise.
with (possibly time-dependent) transition rates 
where 
is the state prior to the jump, while 
is the state after the jump. I’ll assume that this dynamics defines a (continuous-time) Markov process, namely that the numbers 
are the matrix entries of an 
of this process can be characterized by the sequence of occupied states 
and by the times 
at which the transitions 
occur 
:
of each of these transitions to be the logarithmic ratio of transition rates:
by:
is the probability that the system is in state 
. This quantity is also sometimes called the surprisal, as it measures the ‘surprise’ of finding out that the system is in state ![\mathrm{Prob}[\Delta i-\Sigma=-A] = e^{-A}\;\mathrm{Prob}[\Delta i-\Sigma=A]](http://s0.wp.com/latex.php?latex=%5Cmathrm%7BProb%7D%5B%5CDelta+i-%5CSigma%3D-A%5D+%3D+e%5E%7B-A%7D%5C%3B%5Cmathrm%7BProb%7D%5B%5CDelta+i-%5CSigma%3DA%5D&bg=ffffff&fg=333333&s=0 "\mathrm{Prob}[\Delta i-\Sigma=-A] = e^{-A}\;\mathrm{Prob}[\Delta i-\Sigma=A]")
denotes the average over all realizations of the process, then we have the integral fluctuation theorem: 
the ‘entropy production’ or ‘dissipation function’—but how tautological is that?)
and undergoes thermally induced transitions between states 
times the energy difference between the states 
, namely the heat received from the reservoir during the transition. Hence, the mean cumulative skewness 
is nothing but 
with 
), this inequality in turn expresses Landauer’s principle: it impossible to decrease the self-information of the system’s state without dissipating a minimal amount of heat into the environment (in this case 
, the 
in a population, with 
and 
transition rates is completely determined by the 
and 
of 
is a model-dependent function of the effective population size [Sella2005]. Along a given history of mutant fixations, the cumulated skewness 
). This statement may well be the closest thing to a version of the Second Law of Thermodynamics applying to evolutionary dynamics. 
?
be a probability space and 
a measurable involution of 
. Denote 
the pushforward probability measure through this involution, and 
and 
:
by Jensen’s inequality, 
trivially, and 
follows from 
and the change of variables theorem, as follows,
![[0,T]](http://s0.wp.com/latex.php?latex=%5B0%2CT%5D&bg=ffffff&fg=333333&s=0 "[0,T]")
and valued in some (say discrete) state space 
of each state 
satisfying the master equation (aka Kolmogorov equation)
decomposes into ‘variation of self-information’ and ‘cumulative skewness’ along 
:
![\displaystyle{p(\omega)=\pi_{a_{0}}(0)\left[\prod_{j=1}^{N}\phi_{a_{j-1}}(\tau_{j-1},\tau_{j})\gamma_{a_{j-1}a_{j}}(\tau_{j})\right]\phi_{a_{N}}(\tau_{N},T)}](http://s0.wp.com/latex.php?latex=%5Cdisplaystyle%7Bp%28%5Comega%29%3D%5Cpi_%7Ba_%7B0%7D%7D%280%29%5Cleft%5B%5Cprod_%7Bj%3D1%7D%5E%7BN%7D%5Cphi_%7Ba_%7Bj-1%7D%7D%28%5Ctau_%7Bj-1%7D%2C%5Ctau_%7Bj%7D%29%5Cgamma_%7Ba_%7Bj-1%7Da_%7Bj%7D%7D%28%5Ctau_%7Bj%7D%29%5Cright%5D%5Cphi_%7Ba_%7BN%7D%7D%28%5Ctau_%7BN%7D%2CT%29%7D&bg=ffffff&fg=333333&s=0 "\displaystyle{p(\omega)=\pi_{a_{0}}(0)\left[\prod_{j=1}^{N}\phi_{a_{j-1}}(\tau_{j-1},\tau_{j})\gamma_{a_{j-1}a_{j}}(\tau_{j})\right]\phi_{a_{N}}(\tau_{N},T)}")
and 
. It follows from the above lemma that
be defined as above. Then we have
as heat, and Seifert:
