The Azimuth Project
Blog - network theory (part 7)

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

guest post by Jacob Biamonte

This post is part of a series on what John and I like to call Petri net field theory. Stochastic Petri nets can be used to model everything from vending machines to chemical reactions. Chemists have proven some powerful theorems about when these systems have equilibrium states. We’re trying to bind these old ideas into our fancy framework, in hopes that quantum field theory techniques could also be useful in this deep subject. We’ll describe the general theory later; today we’ll do an example from population biology.

Those of you following this series should know that I’m the calculation bunny for this project, with John playing the role of the wolf. If I don’t work quickly, drawing diagrams and trying to keep up with John’s space-bending quasar of information, I’ll be eaten alive! It’s no joke, so please try to respond and pretend to enjoy anything you read here. This will keep me alive for longer. If I did not take notes during our meetings, lots of this stuff would have never made it here, so hope you enjoy.

Amoeba reproduction and competition

Here’s a stochastic Petri net:

It shows a world with one state, amoeba, and two transitions:

reproduction, where one amoeba turns into two. Let’s call the rate constant for this transition α\alpha.

competition, where two amoebas battle for resources and only one survives. Let’s call the rate constant for this transition β\beta.

We are going to analyse this example in several ways. First we’ll study the deterministic dynamics it describes: we’ll look at its rate equation, which turns out to be the logistic equation, familiar in population biology. Then we’ll study the stochastic dynamics, meaning its master equation. That’s where the ideas from quantum field theory come in.

The rate equation

If P(t)P(t) is the population of amoebas at time tt, we can follow the rules explained in Part 3 and crank out this rate equation:

dPdt=αPβP 2 \frac{d P}{d t} = \alpha P - \beta P^2

We can rewrite this as

dPdt=kP(1PQ) \frac{d P }{d t}= k P(1-\frac{P}{Q})

where

Q=αβ,k=α Q = \frac{\alpha}{\beta} , \qquad k = \alpha

What’s the meaning of QQ and kk?

QQ is the carrying capacity, that is, the maximum sustainable population the environment can support.

kk is the growth rate describing the approximately exponential growth of population when P(t)P(t) is small.

It’s a rare treat to find such an important differential equation that can be solved by analytical methods. Let’s enjoy solving it.

We start by separating variables and integrating both sides:

dPP(1P/Q)=kdt\int \frac{d P}{P (1-P/Q)} = \int k d t

We need to use partial fractions on the left side above, resulting in

dPP+dPQP=kdt\int \frac{d P}{P} + \int \frac{d P}{Q-P} = \int k d t

and so we pick up a constant CC, let A=±e CA=\pm e^{-C}, and rearrange things as

QPP=Ae kt \frac{Q-P}{P}=A e^{-k t}

and the population as a function of time becomes

P(t)=Q1+Ae kt P(t) = \frac{Q}{1+A e^{-k t}}

At t=0t=0 we can determine AA uniquely. We write P 0:=P(0)P_0 := P(0) and AA becomes

A=QP 0P 0 A = \frac{Q-P_0}{P_0}

The model now becomes very intuitive. Let’s set Q=k=1Q = k=1 and make a plot for various values of AA:

We arrive at three distinct cases:

equilibrium (A=0A=0). The horizontal blue line corresponds to the case where the initial population P 0P_0 exactly equals the carrying capacity. In this case the population is constant.

dieoff (A<0A \lt 0). The three decaying curves above the horizontal blue line correspond to cases where initial population is higher than the carrying capacity. The population dies off over time and approaches the carrying capacity.

growth (A>0A \gt 0). The four increasing curves below the horizontal blue line represent cases where the initial population is lower than the carrying capacity. Now the population grows over time and approaches the carrying capacity.

The master equation

Next, let us follow the rules explained in Part 6 to write down the master equation for our example. Remember, now we write:

Ψ(t)= n=0ψ n(t)z n \Psi(t) = \sum_{n = 0} \psi_n(t) z^n

where ψ n(t)\psi_n(t) is the probability of having nn amoebas at time tt, and zz is a formal variable. The master equation says

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

where HH is an operator on formal power series called the Hamiltonian. To get the Hamiltonian we take each transition in our Petri net and build an operator built from creation and annihilation operators, as follows. Reproduction works like this:

while competition works like this:

Here aa is the annihilation operator, a a^\dagger is the creation operator and N=a aN = a^\dagger a is the number operator. Last time John explained precisely how the NN‘s arise. So the theory is already in place, and we arrive at this Hamiltonian:

H=α(a a aN)+β(a aaN(N1)) H = \alpha (a^\dagger a^\dagger a - N) \;\; + \;\; \beta(a^\dagger a a - N(N-1))

Remember, α\alpha is the rate constant for reproduction, while β\beta is the rate constant for competition.

The master equation can be solved: it’s equivalent to ddt(e tHΨ(t))=0\frac{d}{d t}(e^{-t H}\Psi(t))=0 so that e tHΨ(t)e^{-t H}\Psi(t) is constant, and so

Ψ(t)=e tHΨ(0) \Psi(t) = e^{t H}\Psi(0)

and that’s it! We can calculate the time evolution starting from any initial probability distribution of populations. Maybe everyone is already used to this, but I find it rather remarkable.

Here’s how it works. We pick a population, say nn amoebas at t=0t=0. This would mean Ψ(0)=z n\Psi(0) = z^n. We then evolve this state using e tHe^{t H}. We will expand this operator as

e tH = n=01n!H nt n = 1+tH+12t 2H n+ \begin{aligned} e^{t H} &=& \sum_{n=0} \frac{1}{n!} H^n t^n \\ \\ &=& 1 + t H + \frac{1}{2}t^2 H^n + \cdots \end{aligned}

This operator contains the full information for the evolution of the system. It contains the histories of all possible amoeba populations—an amoeba mosaic if you will. From this, we can construct amoeba Feynman diagrams.

To do this, we work out each of the H nH^n terms in the expansion above. The first-order terms correspond to the Hamiltonian acting once. These are proportional to either α\alpha or β\beta. The second-order terms correspond to the Hamiltonian acting twice. These are proportional to either α 2\alpha^2, αβ\alpha\beta or β 2\beta^2. And so on.

This is where things start to get interesting! We will consider two possibilities for the second-order terms:

1) A lone amoeba (Ψ(0)=z\Psi(0) = z) reproduces and splits into two. In the battle of the century, the resulting amoebas compete and one dies. At the end we have:

αβ2(a aa)(a a a)z \frac{\alpha \beta}{2} (a^\dagger a a)(a^\dagger a^\dagger a) z

We can draw this as a Feynman diagram:

You might find this tale grim, and you may not like the odds either. It’s true, the odds could be better, but people are worse off than amoebas! The great Japanese swordsman Miyamoto Musashi quoted the survival odds of fair sword duels as 1/3, seeing that 1/3 of the time both participants die. A remedy is to cheat, but these amoeba are competing honestly.

2) We start with two amoebas, so the initial state is Ψ(0)=z 2\Psi(0) = z^2. One of these amoebas splits into two. One of these then gets into an argument with the original amoeba over the Azimuth blog. The amoeba who solved all John’s puzzles survived. At the end we have

αβ2(a aa)(a a a)z 2 \frac{\alpha \beta}{2} (a^\dagger a a)(a^\dagger a^\dagger a) z^2

with corresponding Feynman diagram:

This should give an idea of how this all works. The exponential of the Hamiltonian gives all possible histories, and each of these can be translated into a Feynman diagram. In a future blog entry, we might explain this theory in detail.

An equilibrium state

We’ve seen the equilibrium solution for the rate equation; now let’s look for equilibrium solutions of the master equation. 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 of stochastic Petri nets, there exists an equilibrium solution 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 knowing how many things are in one state tells you nothing about how many are in another!

Here’s a nice quote from this paper:

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.

The ‘deficiency zero theorem’ is a result of Feinberg, which says that for a large class of stochastic Petri nets, the rate equation has a unique equilibrium solution. Anderson, Craciun and Kurtz showed how to use this fact to get equilibrium solutions of the master equation!

We will consider this in future posts. For now, we need to talk a bit about ‘coherent states’.

These are all over the place in quantum theory. Legend (or at least Wikipedia) has it that Erwin Schrödinger himself discovered coherent states when he was looking for states of a quantum system that look ‘as classical as possible’. Suppose you have a quantum harmonic oscillator. Then the uncertainty principle says that

ΔpΔq/2 \Delta p \Delta q \ge \hbar/2

where Δp\Delta p is the uncertainty in the momentum and Δq\Delta q is the uncertainty in position. Suppose we want to make ΔpΔq\Delta p \Delta q as small as possible, and suppose we also want Δp=Δq\Delta p = \Delta q. Then we need our particle to be in a ‘coherent state’. That’s the definition. For the quantum harmonic oscillator, there’s a way to write quantum states as formal power series

Ψ= n=0ψ nz n \Psi = \sum_{n = 0} \psi_n z^n

where ψ n\psi_n is the amplitude for having nn quanta of energy. A coherent state then looks like this:

Ψ=e cz= n=0c nn!z n \Psi = e^{c z} = \sum_{n = 0} \frac{c^n}{n!} z^n

where cc can be any complex number. Here we have omitted a constant factor necessary to normalize the state.

We can also use coherent states in classical stochastic systems like collections of amoebas! Now the coefficient of z nz^n tells us the probability of having nn amoebas, so cc had better be real. And probabilities should sum to 1, so we really should normalize Ψ\Psi as follows:

Ψ=e cze c=e c n=0c nn!z n \Psi = \frac{e^{c z}}{e^c} = e^{-c} \sum_{n = 0} \frac{c^n}{n!} z^n

Now, the probability distribution

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

is called Poisson distribution. So, for starters you can think of a ‘coherent state’ as an over-educated way of talking about a Poisson distribution.

Let’s work out the expected number of amoebas in this Poisson distribution. In the answers to the puzzles in Part 6, we started using this abbreviation

Ψ= n=0 ψ n \sum \Psi = \sum_{n = 0}^\infty \psi_n

We also saw that the expected number of amoebas in the probability distribution Ψ\Psi is

NΨ \sum N \Psi

What does this equal? Remember that N=a aN = a^\dagger a. The annihilation operator aa is just ddz\frac{d}{d z}, so

aΨ=cΨ a \Psi = c \Psi

and we get

NΨ=a aΨ=ca Ψ \sum N \Psi = \sum a^\dagger a \Psi = c \sum a^\dagger \Psi

But we saw in Part 5 that a a^\dagger is stochastic, meaning

a Ψ=Ψ\sum a^\dagger \Psi = \sum \Psi

for any Ψ\Psi. Furthermore, our Ψ\Psi here has

Ψ=1\sum \Psi = 1

since it’s a probability distribution. So:

NΨ=ca Ψ=cΨ=c \sum N \Psi = c \sum a^\dagger \Psi = c \sum \Psi = c

The expected number of amoebas is just cc!

Puzzle 1. This calculation must be wrong if cc is negative: there can’t be a negative number of amoebas. What goes wrong then?

Puzzle 2. Use the same tricks to calculate the standard deviation of the number of amoebas in the Poisson distribution Ψ\Psi.

Now let’s return to our problem and consider the initial amoeba state

Ψ=e cz \Psi = e^{c z}

Here aren’t bothering to normalize it, because we’re going to look for equilibrium solutions to the master equation, meaning solutions where Ψ(t)\Psi(t) doesn’t change with time. So, we want to solve

HΨ=0 H \Psi = 0

Since this equation is linear, the normalization of Ψ\Psi doesn’t really matter.

Remember,

HΨ=α(a a aN)Ψ+β(a aaN(N1))Ψ H\Psi = \alpha (a^\dagger a^\dagger a - N)\Psi + \beta(a^\dagger a a - N(N-1)) \Psi

Let’s work this out. First consider the two α\alpha terms:

a a aΨ=cz 2Ψ a^\dagger a^\dagger a \Psi = c z^2 \Psi

and

NΨ=a aΨ=czΨ -N \Psi = -a^\dagger a\Psi = -c z \Psi

Likewise for the β\beta terms we find

a aaΨ=c 2zΨa^\dagger a a\Psi=c^2 z \Psi

and

N(N1)ψ=a a aaΨ=c 2z 2Ψ-N(N-1)\psi = -a^\dagger a^\dagger a a \Psi = -c^2 z^2\Psi

Here I’m using something John showed in Part 6: the product a 2a 2{a^\dagger}^2 a^2 equals the ‘falling power’ N(N1)N(N-1).

The sum of all four terms must vanish. This happens whenever

α(cz 2cz)+β(c 2zc 2z 2)=0\alpha(c z^2 - c z)+\beta(c^2 z-c^2 z^2) = 0

which is satisfied for

c=αβc= \frac{\alpha}{\beta}

Yipee! We’ve found an equilibrium solution, since we found a value for cc that makes HΨ=0H \Psi = 0. Even better, we’ve seen that the expected number of amoebas in this equilibrium state is

αβ \frac{\alpha}{\beta}

This is just the same as the equilibrium population we saw in the rate equation—that is, the logistic equation! That’s pretty cool, but it’s no coincidence: in fact, Anderson proved it works like this for lots of stochastic Petri nets.

I’m not sure what’s up next or what’s in store, since I’m blogging at gun point from inside a rabbit cage.

Give me all your blog posts, get out of that rabbit cage and reach for the sky!

I’d imagine we’re going to work out the theory behind this example and prove the existence of equilibrium solutions in master equations more generally. One idea John had was to have me start a night shift—that way you’ll get Azimuth posts 24 hours a day.

category: blog