The Azimuth Project
Blog - network theory (part 9) (changes)

Showing changes from revision #13 to #14: Added | Removed | Changed

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

jointly written with Brendan Fong

Last time we reviewed the rate equation and the master equation. Both of them describe reactions where things of various kinds can react and turn into other things. But:

• In the rate equation, we assume the number of things varies continuously and is known precisely.

• In the master equation, we assume the number of things varies discretely and is known only probabilistically.

This should remind you of the difference between classical mechanics and quantum mechanics. However, the master equation is not quantum, it’s stochastic: it involves probabilities, but there’s no uncertainty principle going on. Still, a lot of the math is similar.

Now, given an equilibrium solution to the rate equation—one that doesn’t change with time—we’ll try to find a solution to the master equation with the same property. We won’t always succeed—but we often can! The theorem saying how was proved here:

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

We will translate their proof into the language of annihilation and creation operators. This emphasizes the analogy to quantum mechanics. In particular, our equilibrium solution of the master equation is just like what people call a ‘coherent state’ in quantum mechanics.

So, if you know about quantum mechanics and coherent states, you should be happy. But if you don’t, fear not!—we won’t assume any knowledge of that stuff.

The rate equation

To construct our equilibrium solution of the master equation, we need a special type of solution to our rate equation. We call this type a ‘complex balanced solution’. This essentially means that not only is the net rate of production of each species zero, but the net rate of production of each possible bunch of species of zero.

Before we make this more precise, let’s remind ourselves of the basic setup.

We’ll consider a stochastic Petri net with a finite set SS of species and a finite set TT of transitions. For convenience let’s take S={1,,k}S = \{1,\dots, k\}, so our species are numbered from 1 to kk. Then each transition τ\tau has an input vector m(τ) km(\tau) \in \mathbb{N}^k and output vector n(τ) kn(\tau) \in \mathbb{N}^k. These say how many things of each species go in, and how many go out. Each transition also has rate constant r(τ)[0,)r(\tau) \in [0,\infty), which says how rapidly it happens.

The rate equation concerns a vector x(t)[0,) kx(t) \in [0,\infty)^k whose iith component is the number of things of the iith species at time tt. Note: we’re assuming this number of things varies continuously and is known precisely! This should remind you of classical mechanics. So, we’ll call x(t)x(t), or indeed any vector in [0,) k[0,\infty)^k, a classical state.

The rate equation says how the classical state x(t)x(t) changes with time:

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

You may wonder what x m(τ)x^{m(\tau)} means: after all, we’re taking a vector to a vector power! It’s just an abbreviation, which we’ve seen plenty of times before. If x kx \in \mathbb{R}^k is a list of numbers and m km \in \mathbb{N}^k is a list of natural numbers, we define

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

We’ll also use this notation when xx is a list of operators.

Complex balance

The vectors m(τ)m(\tau) and n(τ)n(\tau) are examples of what chemists call complexes. A complex is just a bunch of things of each species. For example, if the set SS consists of three species, the complex (1,0,5)(1,0,5) can be thought of as a bunch consisting of one thing of the first species, none of the second species, and five of the third species.

For our Petri net, the set of complexes is the set k\mathbb{N}^k, and the complexes of particular interest are the input complex m(τ)m(\tau) and the output complex n(τ)n(\tau) of each transition τ\tau.

We shall say a classical state c[0,) kc \in [0,\infty)^k is complex balanced if for all complexes κ k\kappa \in \mathbb{N}^k we have

(1) {τ:m(τ)=κ}r(τ)c m(τ)= {τ:n(τ)=κ}r(τ)c m(τ) \sum_{\{\tau : m(\tau) = \kappa\}} r(\tau) c^{m(\tau)} =\sum_{\{\tau : n(\tau) = \kappa\}} r(\tau) c^{m(\tau)}

The left hand side of the above equation, which sums over the transitions in which the input complex is κ\kappa, gives the rate of consumption of the complex κ\kappa. The right hand side, which sums over the transitions in which the output complex is κ\kappa, gives the rate of production of κ\kappa. Thus the above equation requires that the net rate of production of the complex κ\kappa is zero when the number of things of each species is given by the vector cc.

Puzzle 1. Show that if a classical state c[0,) kc \in [0,\infty)^k is complex balanced, and we set x(t)=cx(t) = c for all tt, then x(t)x(t) is a solution of the rate equation.

Since x(t)x(t) doesn’t change with time here, we call it an equilibrium solution of the rate equation. Since x(t)=cx(t) = c is complex balanced, we call it complex balanced equilibrium solution.

The master equation

We’ve seen that any complex balanced classical state gives an equilibrium solution of the rate equation. The Anderson–Craciun–Kurtz Theorem says that it also gives an equilibrium solution of the master equation.

The master equation concerns a formal power series

Ψ(t)= n kψ n(t)z n \Psi(t) = \sum_{n \in \mathbb{N}^k} \psi_n(t) z^n

where

z n=z 1 n 1z k n k z^n = z_1^{n_1} \cdots z_k^{n_k}

and

ψ n(t)=ψ n 1,,n k(t) \psi_n(t) = \psi_{n_1, \dots,n_k}(t)

is the probability that at time tt our system has n 1n_1 things of the first species, n 2n_2 of the second species, and so on.

Note: now we’re assuming this number of things varies discretely and is known only probabilistically! So, we’ll call Ψ(t)\Psi(t), or indeed any formal power series of this sort where the coefficients are probabilities summing to 1, a stochastic state. (Earlier we just called it a ‘state’, but that would get confusing now!)

The master equation says how the stochastic state Ψ(t)\Psi(t) changes with time:

ddtΨ(t)=HΨ(t) \frac{d}{d t} \Psi(t) = H \Psi(t)

where the Hamiltonian HH is:

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

The notation here is designed to neatly summarize some big products of annihilation and creation operators. For any vector n kn \in \mathbb{N}^k, we have

a n=a 1 n 1a k n k a^n = a_1^{n_1} \cdots a_k^{n_k}

and

a n=a 1 n 1a k n k {a^\dagger}^n = {a_1^\dagger }^{n_1} \cdots {a_k^\dagger}^{n_k}

Coherent states

Now suppose c[0,) kc \in [0,\infty)^k is a complex balanced equilibrium solution of the rate equation. We want to get an equilibrium solution of the master equation. How do we do it?

For any c[0,) kc \in [0,\infty)^k there is a stochastic state called a coherent state, defined by

Ψ c=e cze c \Psi_c = \frac{e^{c z}}{e^c}

Here we are using some very terse abbreviations. Namely, we are defining

e c=e c 1e c k e^{c} = e^{c_1} \cdots e^{c_k}

and

e cz=e c 1z 1e c kz k e^{c z} = e^{c_1 z_1} \cdots e^{c_k z_k}

Equivalently,

e cz= n kc nn!z n e^{c z} = \sum_{n \in \mathbb{N}^k} \frac{c^n}{n!}z^n

where c nc^n and z nz^n are defined as products in our usual way, and

n!=n 1!n k! n! = n_1! \, \cdots \, n_k!

Either way, if you unravel the abbrevations, here’s what you get:

Ψ c=e (c 1++c k) n kc 1 n 1c k n kn 1!n k!z 1 n 1z k n k \Psi_c = e^{-(c_1 + \cdots + c_k)} \, \sum_{n \in \mathbb{N}^k} \frac{c_1^{n_1} \cdots c_k^{n_k}} {n_1! \, \cdots \, n_k! } \, z_1^{n_1} \cdots z_k^{n_k}

Maybe now you see why we like the abbreviations.

The name ‘coherent state’ comes from quantum mechanics. In quantum mechanics, we think of a coherent state Ψ c\Psi_c as the ‘quantum state’ that best approximates the classical state cc. But we’re not doing quantum mechanics now, we’re doing probability theory. Ψ c\Psi_c isn’t a ‘quantum state’, it’s a stochastic state.

In probability theory, people like Poisson distributions. In the state Ψ c\Psi_c, the probability of having n in_i things of the iith species is equal to

e c ic i n in i! e^{-c_i} \, \frac{c_i^{n_i}}{n_i!}

This is precisely the definition of a Poisson distribution with mean equal to c ic_i. We can multiply a bunch of factors like this, one for each species, to get

e cc nn! e^{-c} \frac{c^n}{n!}

This is the probability of having n 1n_1 things of the first species, n 2n_2 things of the second, and so on, in the state Ψ c\Psi_c. So, the state Ψ c\Psi_c is a product of independent Poisson distributions. In particular, knowing how many things there are of one species says nothing all about how many things there are of any other species!

It is remarkable that such a simple state can give an equilibrium solution of the master equation, even for very complicated stochastic Petri nets. But it’s true—at least if cc is complex balanced.

The Anderson–Craciun–Kurtz Theorem

Now we’re ready to state and prove the big result:

Theorem (Anderson–Craciun–Kurtz). Suppose c[0,) kc \in [0,\infty)^k is a complex balanced equilibrium solution of the rate equation. Then HΨ c=0H \Psi_c = 0.

It follows that Ψ c\Psi_c is an equilibrium solution of the master equation. In other words, if we take Ψ(t)=Ψ c\Psi(t) = \Psi_c for all times tt, the master equation holds:

ddtΨ(t)=HΨ(t) \frac{d}{d t} \Psi(t) = H \Psi(t)

since both sides are zero.

Proof. To prove the Anderson–Craciun–Kurtz Theorem, we just need to show that HΨ c=0H \Psi_c = 0. Since Ψ c\Psi_c is a constant times e cze^{c z}, it suffices to show He cz=0H e^{c z} = 0. Remember that

He cz= τTr(τ)(a n(τ)a m(τ))a m(τ)e cz H e^{c z} = \sum_{\tau \in T} r(\tau) \left( {a^\dagger}^{n(\tau)} -{a^\dagger}^{m(\tau)} \right) \, a^{m(\tau)} \, e^{c z}

Since the annihilation operator a ia_i is given by differentiation with respect to z iz_i, while the creation operator a i a^\dagger_i is just multiplying by z iz_i, we have:

He cz= τTr(τ)c m(τ)(z n(τ)z m(τ))e cz H e^{c z} = \sum_{\tau \in T} r(\tau) \, c^{m(\tau)} \left( z^{n(\tau)} - z^{m(\tau)} \right) e^{c z}

Expanding out e cze^{c z} we get:

He cz= i k τTr(τ)c m(τ)(z n(τ)c ii!z iz m(τ)c ii!z i) H e^{c z} = \sum_{i \in \mathbb{N}^k} \sum_{\tau \in T} r(\tau)c^{m(\tau)}\left(z^{n(\tau)}\frac{c^i}{i!}z^i - z^{m(\tau)}\frac{c^i}{i!}z^i\right)

Shifting indices and defining negative powers to be zero:

He cz= i k τTr(τ)c m(τ)(c in(τ)(in(τ))!z ic im(τ)(im(τ))!z i) H e^{c z} = \sum_{i \in \mathbb{N}^k} \sum_{\tau \in T} r(\tau)c^{m(\tau)}\left(\frac{c^{i-n(\tau)}}{(i-n(\tau))!}z^i - \frac{c^{i-m(\tau)}}{(i-m(\tau))!}z^i\right)

So, to show He cz=0H e^{c z} = 0, we need to show this:

i k τTr(τ)c m(τ)c in(τ)(in(τ))!z i= i k τTr(τ)c m(τ)c im(τ)(im(τ))!z i \sum_{i \in \mathbb{N}^k} \sum_{\tau \in T} r(\tau)c^{m(\tau)}\frac{c^{i-n(\tau)}}{(i-n(\tau))!}z^i =\sum_{i \in \mathbb{N}^k} \sum_{\tau \in T} r(\tau)c^{m(\tau)}\frac{c^{i-m(\tau)}}{(i-m(\tau))!}z^i

We do this by splitting the sum over TT according to output and then input complexes, making use of the complex balanced condition:

i k κ k {τ:n(τ)=κ}r(τ)c m(τ)c in(τ)(in(τ))!z i = i k κ kc iκ(iκ)!z i {τ:n(τ)=κ}r(τ)c m(τ) = i k κ kc iκ(iκ)!z i {τ:m(τ)=κ}r(τ)c m(τ) = i k κ k {τ:m(τ)=κ}r(τ)c m(τ)c im(τ)(im(τ))!z i \begin{array}{ccl} \sum_{i \in \mathbb{N}^k} \sum_{\kappa \in \mathbb{N}^k} \sum_{\{\tau : n(\tau)=\kappa\}} r(\tau)c^{m(\tau)}\frac{c^{i-n(\tau)}}{(i-n(\tau))!} \, z^i \; &=& \sum_{i \in \mathbb{N}^k} \sum_{\kappa \in \mathbb{N}^k} \frac{c^{i-\kappa}}{(i-\kappa)!}z^i \sum_{\{\tau : n(\tau) = \kappa\}} r(\tau)c^{m(\tau)} \\ &=& \sum_{i \in \mathbb{N}^k} \sum_{\kappa \in \mathbb{N}^k} \frac{c^{i-\kappa}}{(i-\kappa)!}z^i \sum_{\{\tau : m(\tau) = \kappa\}} r(\tau)c^{m(\tau)} \\ &=& \sum_{i \in \mathbb{N}^k} \sum_{\kappa \in \mathbb{N}^k} \sum_{\{\tau : m(\tau) = \kappa\}} r(\tau)c^{m(\tau)}\frac{c^{i-m(\tau)}}{(i-m(\tau))!} \, z^i \end{array}

This completes the proof!

An example

We’ve already seen one example of this theorem, back in Part 7. We had this stochastic Petri net:

We saw that the rate equation is just the logistic equation, familiar from population biology. The equilibrium solution is complex balanced, because pairs of amoebas are getting created at the same rate as they’re getting destroyed, and single amoebas are getting created at the same rate as they’re getting destroyed. So, the Anderson–Craciun–Kurtz Theorem guarantees that there’s an equilibrium solution of the master equation where the number of amoebas is distributed according to a Poisson distribution. And, we actually found this equilibrium solution!

Next time we’ll look at another example.

category: blog