Network Theory (Part 5)

11 April, 2011

Last time we saw clues that stochastic Petri nets are a lot like quantum field theory, but with probabilities replacing amplitudes. There’s a powerful analogy at work here, which can help us a lot. So, this time I want to make that analogy precise.

But first, let me quickly sketch why it could be worthwhile.

A Poisson process

Consider this stochastic Petri net with rate constant r:

It describes an inexhaustible supply of fish swimming down a river, and getting caught when they run into a fisherman’s net. In any short time \Delta t there’s a chance of about r \Delta t of a fish getting caught. There’s also a chance of two or more fish getting caught, but this becomes negligible by comparison as \Delta t \to 0. Moreover, the chance of a fish getting caught during this interval of time is independent of what happens before or afterwards. This sort of process is called a Poisson process.

Problem. Suppose we start out knowing for sure there are no fish in the fisherman’s net. What’s the probability that he has caught n fish at time t?

At any time there will be some probability of having caught n fish; let’s call this probability \psi(n,t), or \psi(n) for short. We can summarize all these probabilities in a single power series, called a generating function:

\Psi(t) = \sum_{n=0}^\infty \psi(n,t) \, z^n

Here z is a formal variable—don’t ask what it means, for now it’s just a trick. In quantum theory we use this trick when talking about collections of photons rather than fish, but then the numbers \psi(n,t) are complex ‘amplitudes’. Now they are real probabilities, but we can still copy what the physicists do, and use this trick to rewrite the master equation as follows:

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

This describes how the probability of having caught any given number of fish changes with time.

What’s the operator H? Well, in quantum theory we describe the creation of photons using a certain operator on power series called the creation operator:

a^\dagger \Psi = z \Psi

We can try to apply this to our fish. If at some time we’re 100% sure we have n fish, we have

\Psi = z^n

so applying the creation operator gives

a^\dagger \Psi = z^{n+1}

One more fish! That’s good. So, an obvious wild guess is

H = r a^\dagger

where r is the rate at which we’re catching fish. Let’s see how well this guess works.

If you know how to exponentiate operators, you know to solve this equation:

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

It’s easy:

\Psi(t) = \mathrm{exp}(t H) \Psi(0)

Since we start out knowing there are no fish in the net, we have

\Psi(0) = 1

so with our guess for H we get

\Psi(t) = \mathrm{exp}(r t a^\dagger) 1

But a^\dagger is the operator of multiplication by z, so \mathrm{exp}(r t a^\dagger) is multiplication by e^{r t z}, so

\Psi(t) = e^{r t z} = \sum_{n = 0}^\infty \frac{(r t)^n}{n!} \, z^n

So, if our guess is right, the probability of having caught n fish at time t is

\frac{(r t)^n}{n!}

Unfortunately, this can’t be right, because these probabilities don’t sum to 1! Instead their sum is

\sum_{n=0}^\infty \frac{(r t)^n}{n!} = e^{r t}

We can try to wriggle out of the mess we’re in by dividing our answer by this fudge factor. It sounds like a desperate measure, but we’ve got to try something!

This amounts to guessing that the probability of having caught n fish by time t is

\frac{(r t)^n}{n!} \, e^{-r t }

And this is right! This is called the Poisson distribution: it’s famous for being precisely the answer to the problem we’re facing.

So on the one hand our wild guess about H was wrong, but on the other hand it was not so far off. We can fix it as follows:

H = r (a^\dagger - 1)

The extra -1 gives us the fudge factor we need.

So, a wild guess corrected by an ad hoc procedure seems to have worked! But what’s really going on?

What’s really going on is that a^\dagger, or any multiple of this, is not a legitimate Hamiltonian for a master equation: if we define a time evolution operator \exp(t H) using a Hamiltonian like this, probabilities won’t sum to 1! But a^\dagger - 1 is okay. So, we need to think about which Hamiltonians are okay.

In quantum theory, self-adjoint Hamiltonians are okay. But in probability theory, we need some other kind of Hamiltonian. Let’s figure it out.

Probability versus quantum theory

Suppose we have a system of any kind: physical, chemical, biological, economic, whatever. The system can be in different states. In the simplest sort of model, we say there’s some set X of states, and say that at any moment in time the system is definitely in one of these states. But I want to compare two other options:

• In a probabilistic model, we may instead say that the system has a probability \psi(x) of being in any state x \in X. These probabilities are nonnegative real numbers with

\sum_{x \in X} \psi(x) = 1

• In a quantum model, we may instead say that the system has an amplitude \psi(x) of being in any state x \in X. These amplitudes are complex numbers with

\sum_{x \in X} | \psi(x) |^2 = 1

Probabilities and amplitudes are similar yet strangely different. Of course given an amplitude we can get a probability by taking its absolute value and squaring it. This is a vital bridge from quantum theory to probability theory. Today, however, I don’t want to focus on the bridges, but rather the parallels between these theories.

We often want to replace the sums above by integrals. For that we need to replace our set X by a measure space, which is a set equipped with enough structure that you can integrate real or complex functions defined on it. Well, at least you can integrate so-called ‘integrable’ functions—but I’ll neglect all issues of analytical rigor here. Then:

• In a probabilistic model, the system has a probability distribution \psi : X \to \mathbb{R}, which obeys \psi \ge 0 and

\int_X \psi(x) \, d x = 1

• In a quantum model, the system has a wavefunction \psi : X \to \mathbb{C}, which obeys

\int_X | \psi(x) |^2 \, d x= 1

In probability theory, we integrate \psi over a set S \subset X to find out the probability that our systems state is in this set. In quantum theory we integrate |\psi|^2 over the set to answer the same question.

We don’t need to think about sums over sets and integrals over measure spaces separately: there’s a way to make any set X into a measure space such that by definition,

\int_X \psi(x) \, dx = \sum_{x \in X} \psi(x)

In short, integrals are more general than sums! So, I’ll mainly talk about integrals, until the very end.

In probability theory, we want our probability distributions to be vectors in some vector space. Ditto for wave functions in quantum theory! So, we make up some vector spaces:

• In probability theory, the probability distribution \psi is a vector in the space

L^1(X) = \{ \psi: X \to \mathbb{C} \; : \; \int_X |\psi(x)| \, d x < \infty \}

• In quantum theory, the wavefunction \psi is a vector in the space

L^2(X) = \{ \psi: X \to \mathbb{C} \; : \; \int_X |\psi(x)|^2 \, d x < \infty \}

You may wonder why I defined L^1(X) to consist of complex functions when probability distributions are real. I’m just struggling to make the analogy seem as strong as possible. In fact probability distributions are not just real but nonnegative. We need to say this somewhere… but we can, if we like, start by saying they’re complex-valued functions, but then whisper that they must in fact be nonnegative (and thus real). It’s not the most elegant solution, but that’s what I’ll do for now.

Now:

• The main thing we can do with elements of L^1(X), besides what we can do with vectors in any vector space, is integrate one. This gives a linear map:

\int : L^1(X) \to \mathbb{C}

• The main thing we can with elements of L^2(X), besides the besides the things we can do with vectors in any vector space, is take the inner product of two:

\langle \psi, \phi \rangle = \int_X \overline{\psi}(x) \phi(x) \, d x

This gives a map that’s linear in one slot and conjugate-linear in the other:

\langle - , - \rangle :  L^2(X) \times L^2(X) \to \mathbb{C}

First came probability theory with L^1(X); then came quantum theory with L^2(X). Naive extrapolation would say it’s about time for someone to invent an even more bizarre theory of reality based on L^3(X). In this, you’d have to integrate the product of three wavefunctions to get a number! The math of Lp spaces is already well-developed, so give it a try if you want. I’ll stick to L^1 and L^2 today.

Stochastic versus unitary operators

Now let’s think about time evolution:

• In probability theory, the passage of time is described by a map sending probability distributions to probability distributions. This is described using a stochastic operator

U : L^1(X) \to L^1(X)

meaning a linear operator such that

\int U \psi = \int \psi

and

\psi \ge 0 \quad \implies \quad U \psi \ge 0

• In quantum theory the passage of time is described by a map sending wavefunction to wavefunctions. This is described using an isometry

U : L^2(X) \to L^2(X)

meaning a linear operator such that

\langle U \psi , U \phi \rangle = \langle \psi , \phi \rangle

In quantum theory we usually want time evolution to be reversible, so we focus on isometries that have inverses: these are called unitary operators. In probability theory we often consider stochastic operators that are not invertible.

Infinitesimal stochastic versus self-adjoint operators

Sometimes it’s nice to think of time coming in discrete steps. But in theories where we treat time as continuous, to describe time evolution we usually need to solve a differential equation. This is true in both probability theory and quantum theory:

• In probability theory we often describe time evolution using a differential equation called the master equation:

\frac{d}{d t} \psi(t) = H \psi(t)

whose solution is

\psi(t) = \exp(t H)\psi(0)

• In quantum theory we often describe time evolution using a differential equation called Schrödinger’s equation:

i \frac{d}{d t} \psi(t) = H \psi(t)

whose solution is

\psi(t) = \exp(-i t H)\psi(0)

In fact the appearance of i in the quantum case is purely conventional; we could drop it to make the analogy better, but then we’d have to work with ‘skew-adjoint’ operators instead of self-adjoint ones in what follows.

Let’s guess what properties an operator H should have to make \exp(-i t H) unitary for all t. We start by assuming it’s an isometry:

\langle \exp(-i t H) \psi, \exp(-i t H) \phi \rangle = \langle \psi, \phi \rangle

Then we differentiate this with respect to t and set t = 0, getting

\langle -iH \psi, \phi \rangle +  \langle \psi, -iH \phi \rangle = 0

or in other words

\langle H \psi, \phi \rangle = \langle \psi, H \phi \rangle

Physicists call an operator obeying this condition self-adjoint. Mathematicians know there’s more to it, but today is not the day to discuss such subtleties, intriguing though they be. All that matters now is that there is, indeed, a correspondence between self-adjoint operators and well-behaved ‘one-parameter unitary groups’ \exp(-i t H). This is called Stone’s Theorem.

But now let’s copy this argument to guess what properties an operator H must have to make \exp(t H) stochastic. We start by assuming \exp(t H) is stochastic, so

\int \exp(t H) \psi = \int \psi

and

\psi \ge 0 \quad \implies \quad \exp(t H) \psi \ge 0

We can differentiate the first equation with respect to t and set t = 0, getting

\int H \psi = 0

for all \psi.

But what about the second condition,

\psi \ge 0 \quad \implies \quad \exp(t H) \psi \ge 0?

It seems easier to deal with this in the special case when integrals over X reduce to sums. So let’s suppose that happens… and let’s start by seeing what the first condition says in this case.

In this case, L^1(X) has a basis of ‘Kronecker delta functions’: The Kronecker delta function \delta_i vanishes everywhere except at one point i \in X, where it equals 1. Using this basis, we can write any operator on L^1(X) as a matrix.

As a warmup, let’s see what it means for an operator

U: L^1(X) \to L^1(X)

to be stochastic in this case. We’ll take the conditions

\int U \psi = \int \psi

and

\psi \ge 0 \quad \implies \quad U \psi \ge 0

and rewrite them using matrices. For both, it’s enough to consider the case where \psi is a Kronecker delta, say \delta_j.

In these terms, the first condition says

\sum_{i \in X} U_{i j} = 1

for each column j. The second says

U_{i j} \ge 0

for all i, j. So, in this case, a stochastic operator is just a square matrix where each column sums to 1 and all the entries are nonnegative. (Such matrices are often called left stochastic.)

Next, let’s see what we need for an operator H to have the property that \exp(t H) is stochastic for all t \ge 0. It’s enough to assume t is very small, which lets us use the approximation

\exp(t H) = 1 + t H + \cdots

and work to first order in t. Saying that each column of this matrix sums to 1 then amounts to

\sum_{i \in X} \delta_{i j} + t H_{i j} + \cdots  = 1

which requires

\sum_{i \in X} H_{i j} = 0

Saying that each entry is nonnegative amounts to

\delta_{i j} + t H_{i j} + \cdots  \ge 0

When i = j this will be automatic when t is small enough, so the meat of this condition is

H_{i j} \ge 0 \; \mathrm{if} \; i \ne j

So, let’s say H is an infinitesimal stochastic matrix if its columns sum to zero and its off-diagonal entries are nonnegative.

I don’t love this terminology: do you know a better one? There should be some standard term. People here say they’ve seen the term ‘stochastic Hamiltonian’. The idea behind my term is that any infintesimal stochastic operator should be the infinitesimal generator of a stochastic process.

In other words, when we get the details straightened out, any 1-parameter family of stochastic operators

U(t) : L^1(X) \to L^1(X) \qquad   t \ge 0

obeying

U(0) = I

U(t) U(s) = U(t+s)

and continuity:

t_i \to t \quad \implies \quad U(t_i) \psi \to U(t)\psi

should be of the form

U(t) = \exp(t H)

for a unique ‘infinitesimal stochastic operator’ H.

When X is a finite set, this is true—and an infinitesimal stochastic operator is just a square matrix whose columns sum to zero and whose off-diagonal entries are nonnegative. But do you know a theorem characterizing infinitesimal stochastic operators for general measure spaces X? Someone must have worked it out.

Luckily, for our work on stochastic Petri nets, we only need to understand the case where X is a countable set and our integrals are really just sums. This should be almost like the case where X is a finite set—but we’ll need to take care that all our sums converge.

The moral

Now we can see why a Hamiltonian like a^\dagger is no good, while a^\dagger - 1 is good. (I’ll ignore the rate constant r since it’s irrelevant here.) The first one is not infinitesimal stochastic, while the second one is!

In this example, our set of states is the natural numbers:

X = \mathbb{N}

The probability distribution

\psi : \mathbb{N} \to \mathbb{C}

tells us the probability of having caught any specific number of fish.

The creation operator is not infinitesimal stochastic: in fact, it’s stochastic! Why? Well, when we apply the creation operator, what was the probability of having n fish now becomes the probability of having n+1 fish. So, the probabilities remain nonnegative, and their sum over all n is unchanged. Those two conditions are all we need for a stochastic operator.

Using our fancy abstract notation, these conditions say:

\int a^\dagger \psi = \int \psi

and

\psi \ge 0 \; \implies \; a^\dagger \psi \ge 0

So, precisely by virtue of being stochastic, the creation operator fails to be infinitesimal stochastic:

\int a^\dagger \psi \ne 0

Thus it’s a bad Hamiltonian for our stochastic Petri net.

On the other hand, a^\dagger - 1 is infinitesimal stochastic. Its off-diagonal entries are the same as those of a^\dagger, so they’re nonnegative. Moreover:

\int (a^\dagger - 1) \psi = 0

precisely because

\int a^\dagger \psi = \int \psi

You may be thinking: all this fancy math just to understand a single stochastic Petri net, the simplest one of all!

But next time I’ll explain a general recipe which will let you write down the Hamiltonian for any stochastic Petri net. The lessons we’ve learned today will make this much easier. And pondering the analogy between probability theory and quantum theory will also be good for our bigger project of unifying the applications of network diagrams to dozens of different subjects.


Chemistry Puzzle

7 April, 2011

Three elements were named after mischievous sprites. Two of them are real, while the third was just a mischievous trick played by Mother Nature herself. What are these elements, how did they get their names, and what was the trick?


Network Theory (Part 4)

6 April, 2011

Last time I explained the rate equation of a stochastic Petri net. But now let’s get serious: let’s see what’s really stochastic— that is, random—about a stochastic Petri net. For this we need to forget the rate equation (temporarily) and learn about the ‘master equation’. This is where ideas from quantum field theory start showing up!

A Petri net has a bunch of states and a bunch of transitions. Here’s an example we’ve already seen, from chemistry:

The states are in yellow, the transitions in blue. A labelling of our Petri net is a way of putting some number of things in each state. We can draw these things as little black dots:

In this example there are only 0 or 1 things in each state: we’ve got one atom of carbon, one molecule of oxygen, one molecule of sodium hydroxide, one molecule of hydrochloric acid, and nothing else. But in general, we can have any natural number of things in each state.

In a stochastic Petri net, the transitions occur randomly as time passes. For example, as time passes we could see a sequence of transitions like this:

 

 

 

 

 

 

Each time a transition occurs, the number of things in each state changes in an obvious way.

The Master Equation

Now, I said the transitions occur ‘randomly’, but that doesn’t mean there’s no rhyme or reason to them! The miracle of probability theory is that it lets us state precise laws about random events. The law governing the random behavior of a stochastic Petri net is called the ‘master equation’.

In a stochastic Petri net, each transition has a rate constant, a nonnegative real number. Roughly speaking, this determines the probability of that transition.

A bit more precisely: suppose we have a Petri net that is labelled in some way at some moment. Then the probability that a given transition occurs in a short time \Delta t is approximately:

• the rate constant for that transition, times

• the time \Delta t, times

• the number of ways the transition can occur.

More precisely still: this formula is correct up to terms of order (\Delta t)^2. So, taking the limit as \Delta t \to 0, we get a differential equation describing precisely how the probability of the Petri net having a given labelling changes with time! And this is the master equation.

Now, you might be impatient to actually see the master equation, but that would be rash. The true master doesn’t need to see the master equation. It sounds like a Zen proverb, but it’s true. The raw beginner in mathematics wants to see the solutions of an equation. The more advanced student is content to prove that the solution exists. But the master is content to prove that the equation exists.

A bit more seriously: what matters is understanding the rules that inevitably lead to some equation: actually writing it down is then straightforward.

And you see, there’s something I haven’t explained yet: “the number of ways the transition can occur”. This involves a bit of counting. Consider, for example, this Petri net:

Suppose there are 10 rabbits and 5 wolves.

• How many ways can the birth transition occur? Since birth takes one rabbit as input, it can occur in 10 ways.

• How many ways can predation occur? Since predation takes one rabbit and one wolf as inputs, it can occur in 10 × 5 = 50 ways.

• How many ways can death occur? Since death takes one wolf as input, it can occur in 5 ways.

Or consider this one:

Suppose there are 10 hydrogen atoms and 5 oxygen atoms. How many ways can they form a water molecule? There are 10 ways to pick the first hydrogen, 9 ways to pick the second hydrogen, and 5 ways to pick the oxygen. So, there are

10 \times 9 \times 5 = 450

ways.

Note that we’re treating the hydrogen atoms as distinguishable, so there are 10 \times 9 ways to pick them, not \frac{10 \times 9}{2} = \binom{9}{2}. In general, the number of ways to choose M distinguishable things from a collection of L is the falling power

L^{\underline{M}} = L \cdot (L - 1) \cdots (L - M + 1)

where there are M factors in the product, but each is 1 less than the preceding one—hence the term ‘falling’.

Okay, now I’ve given you all the raw ingredients to work out the master equation for any stochastic Petri net. The previous paragraph was a big fat hint. One more nudge and you’re on your own:

Puzzle. Suppose we have a stochastic Petri net with k states and just one transition, whose rate constant is r. Suppose the ith state appears m_i times as the input of this transition and n_i times as the output. A labelling of this stochastic Petri net is a k-tuple of natural numbers \ell = (\ell_1, \dots, \ell_k) saying how many things are in each state. Let \psi_\ell(t) be the probability that the labelling is \ell at time t. Then the master equation looks like this:

\frac{d}{d t} \psi_{\ell}(t)  = \sum_{\ell'} H_{\ell \ell'} \psi_{\ell'}(t)

for some matrix of real numbers H_{\ell \ell'}. What is this matrix?

You can write down a formula for this matrix using what I’ve told you. And then, if you have a stochastic Petri net with more transitions, you can just compute the matrix for each transition using this formula, and add them all up.

Someday I’ll tell you the answer to this puzzle, but I want to get there by a strange route: I want to guess the master equation using ideas from quantum field theory!

Some clues

Why?

Well, if we think about a stochastic Petri net whose labelling undergoes random transitions as I’ve described, you’ll see that any possible ‘history’ for the labelling can be drawn in a way that looks like a Feynman diagram. In quantum field theory, Feynman diagrams show how things interact and turn into other things. But that’s what stochastic Petri nets do, too!

For example, if our Petri net looks like this:

then a typical history can be drawn like this:

Some rabbits and wolves come in on top. They undergo some transitions as time passes, and go out on the bottom. The vertical coordinate is time, while the horizontal coordinate doesn’t really mean anything: it just makes the diagram easier to draw.

If we ignore all the artistry that makes it cute, this Feynman diagram is just a graph with states as edges and transitions as vertices. Each transition occurs at a specific time.

We can use these Feynman diagrams to compute the probability that if we start it off with some labelling at time t_1, our stochastic Petri net will wind up with some other labelling at time t_2. To do this, we just take a sum over Feynman diagrams that start and end with the given labellings. For each Feynman diagram, we integrate over all possible times at which the transitions occur. And what do we integrate? Just the product of the rate constants for those transitions!

That was a bit of a mouthful, and it doesn’t really matter if you followed it in detail. What matters is that it sounds a lot like stuff you learn when you study quantum field theory!

That’s one clue that something cool is going on here. Another is the master equation itself:

\frac{d}{dt} \psi_{\ell}(t) = \sum_{\ell'} H_{\ell \ell'} \psi_{\ell'}(t)

This looks a lot like Schrödinger’s equation, the basic equation describing how a quantum system changes with the passage of time.

We can make it look even more like Schrödinger’s equation if we create a vector space with the labellings \ell as a basis. The numbers \psi_\ell(t) will be the components of some vector \psi(t) in this vector space. The numbers H_{\ell \ell'} will be the matrix entries of some operator H on that vector space. And the master equation becomes:

\frac{d}{d t} \psi(t) = H \psi(t)

Compare Schrödinger’s equation:

i \frac{d}{d t} \psi(t) = H \psi(t)

The only visible difference is that factor of i!

But of course this is linked to another big difference: in the master equation \psi describes probabilities, so it’s a vector in a real vector space. In quantum theory \psi describes probability amplitudes, so it’s a vector in a complex Hilbert space.

Apart from this huge difference, everything is a lot like quantum field theory. In particular, our vector space is a lot like the Fock space one sees in quantum field theory. Suppose we have a quantum particle that can be in k different states. Then its Fock space is the Hilbert space we use to describe an arbitrary collection of such particles. It has an orthonormal basis denoted

| \ell_1 \cdots \ell_k \rangle

where \ell_1, \dots, \ell_k are natural numbers saying how many particles there are in each state. So, any vector in Fock space looks like this:

\psi = \sum_{\ell_1, \dots, \ell_k} \; \psi_{\ell_1 , \dots, \ell_k} \,  | \ell_1 \cdots \ell_k \rangle

But if write the whole list \ell_1, \dots, \ell_k simply as \ell, this becomes

\psi = \sum_{\ell} \psi_{\ell}   | \ell \rangle

This is almost like what we’ve been doing with Petri nets!—except I hadn’t gotten around to giving names to the basis vectors.

In quantum field theory class, I learned lots of interesting operators on Fock space: annihilation and creation operators, number operators, and so on. So, when I bumped into this master equation

\frac{d}{d t} \psi(t) = H \psi(t)

it seemed natural to take the operator H and write it in terms of these. There was an obvious first guess, which didn’t quite work… but thinking a bit harder eventually led to the right answer. Later, it turned out people had already thought about similar things. So, I want to explain this.

When I first started working on this stuff, I was focused on the difference between collections of indistinguishable things, like bosons or fermions, and collections of distinguishable things, like rabbits or wolves.

But with the benefit of hindsight, it’s even more important to think about the difference between ordinary quantum theory, which is all about probability amplitudes, and the game we’re playing now, which is all about probabilities. So, next time I’ll explain how we need to modify quantum theory so that it’s about probabilities. This will make it easier to guess a nice formula for H.


Network Theory (Part 3)

3 April, 2011

As we saw last time, a Petri net is a picture that shows different kinds of things and processes that turn bunches of things into other bunches of things, like this:

The kinds of things are called states and the processes are called transitions. We see such transitions in chemistry:

H + OH → H2O

and population biology:

amoeba → amoeba + amoeba

and the study of infectious diseases:

infected + susceptible → infected + infected

and many other situations.

A “stochastic” Petri net says the rate at which each transition occurs. We can think of these transitions as occurring randomly at a certain rate—and then we get a stochastic process described by something called the “master equation”. But for starters, we’ve been thinking about the limit where there are very many things in each state. Then the randomness washes out, and the expected number of things in each state changes deterministically in a manner described by the “rate equation”.

It’s time to explain the general recipe for getting this rate equation! It looks complicated at first glance, so I’ll briefly state it, then illustrate it with tons of examples, and then state it again.

One nice thing about stochastic Petri nets is that they let you dabble in many sciences. Last time we got a tiny taste of how they show up in population biology. This time we’ll look at chemistry and models of infectious diseases. I won’t dig very deep, but take my word for it: you can do a lot with stochastic Petri nets in these subjects! I’ll give some references in case you want to learn more.

Rate equations: the general recipe

Here’s the recipe, really quick:

A stochastic Petri net has a set of states and a set of transitions. Let’s concentrate our attention on a particular transition. Then the ith state will appear m_i times as the input to that transition, and n_i times as the output. Our transition also has a reaction rate 0 \le r < \infty.

The rate equation answers this question:

\frac{d x_i}{d t} = ???

where x_i(t) is the number of things in the ith state at time t. The answer is a sum of terms, one for each transition. Each term works the same way. For the transition we’re looking at, it’s

r (n_i - m_i) x_1^{m_1} \cdots x_k^{m_k}

The factor of (n_i - m_i) shows up because our transition destroys m_i things in the ith state and creates n_i of them. The big product over all states, x_1^{m_1} \cdots x_k^{m_k} , shows up because our transition occurs at a rate proportional to the product of the numbers of things it takes as inputs. The constant of proportionality is the reaction rate r.

The formation of water (1)

But let’s do an example. Here’s a naive model for the formation of water from atomic hydrogen and oxygen:

This Petri net has just one transition: two hydrogen atoms and an oxygen atom collide simultaneously and form a molecule of water. That’s not really how it goes… but what if it were? Let’s use [\mathrm{H}] for the number of hydrogen atoms, and so on, and let the reaction rate be \alpha. Then we get this rate equation:

\begin{array}{ccl}  \frac{d [\mathrm{H}]}{d t} &=& - 2 \alpha [\mathrm{H}]^2 [\mathrm{O}]   \\  \\  \frac{d [\mathrm{O}]}{d t} &=& - \alpha [\mathrm{H}]^2 [\mathrm{O}]   \\  \\  \frac{d [\mathrm{H}_2\mathrm{O}]}{d t} &=& \alpha [\mathrm{H}]^2 [\mathrm{O}] \end{array}

See how it works? The reaction occurs at a rate proportional to the product of the numbers of things that appear as inputs: two H’s and one O. The constant of proportionality is the rate constant \alpha. So, the reaction occurs at a rate equal to \alpha [\mathrm{H}]^2 [\mathrm{O}] . Then:

• Since two hydrogen atoms get used up in this reaction, we get a factor of -2 in the first equation.

• Since one oxygen atom gets used up, we get a factor of -1 in the second equation.

• Since one water molecule is formed, we get a factor of +1 in the third equation.

The formation of water (2)

Let me do another example, just so chemists don’t think I’m an absolute ninny. Chemical reactions rarely proceed by having three things collide simultaneously—it’s too unlikely. So, for the formation of water from atomic hydrogen and oxygen, there will typically be an intermediate step. Maybe something like this:

Here OH is called a ‘hydroxyl radical’. I’m not sure this is the most likely pathway, but never mind—it’s a good excuse to work out another rate equation. If the first reaction has rate constant \alpha and the second has rate constant \beta, here’s what we get:

\begin{array}{ccl}  \frac{d [\mathrm{H}]}{d t} &=& - \alpha [\mathrm{H}] [\mathrm{O}] - \beta [\mathrm{H}] [\mathrm{OH}]   \\  \\  \frac{d [\mathrm{OH}]}{d t} &=&  \alpha [\mathrm{H}] [\mathrm{O}] -  \beta [\mathrm{H}] [\mathrm{OH}]   \\  \\  \frac{d [\mathrm{O}]}{d t} &=& - \alpha [\mathrm{H}] [\mathrm{O}]   \\  \\  \frac{d [\mathrm{H}_2\mathrm{O}]}{d t} &=&   \beta [\mathrm{H}] [\mathrm{OH}]  \end{array}

See how it works? Each reaction occurs at a rate proportional to the product of the numbers of things that appear as inputs. We get minus signs when a reaction destroys one thing of a given kind, and plus signs when it creates one. We don’t get factors of 2 as we did last time, because now no reaction creates or destroys two of anything.

The dissociation of water (1)

In chemistry every reaction comes with a reverse reaction. So, if hydrogen and oxygen atoms can combine to form water, a water molecule can also ‘dissociate’ into hydrogen and oxygen atoms. The rate constants for the reverse reaction can be different than for the original reaction… and all these rate constants depend on the temperature. At room temperature, the rate constant for hydrogen and oxygen to form water is a lot higher than the rate constant for the reverse reaction. That’s why we see a lot of water, and not many lone hydrogen or oxygen atoms. But at sufficiently high temperatures, the rate constants change, and water molecules become more eager to dissociate.

Calculating these rate constants is a big subject. I’m just starting to read this book, which looked like the easiest one on the library shelf:

• S. R. Logan, Chemical Reaction Kinetics, Longman, Essex, 1996.

But let’s not delve into these mysteries today. Let’s just take our naive Petri net for the formation of water and turn around all the arrows, to get the reverse reaction:

If the reaction rate is \alpha, here’s the rate equation:

\begin{array}{ccl}  \frac{d [\mathrm{H}]}{d t} &=& 2 \alpha [\mathrm{H}_2\mathrm{O}]   \\  \\  \frac{d [\mathrm{O}]}{d t} &=& \alpha [\mathrm{H}_2 \mathrm{O}]   \\  \\  \frac{d [\mathrm{H}_2\mathrm{O}]}{d t} &=& - \alpha [\mathrm{H}_2 \mathrm{O}]   \end{array}

See how it works? The reaction occurs at a rate proportional to [\mathrm{H}_2\mathrm{O}], since it has just a single water molecule as input. That’s where the \alpha [\mathrm{H}_2\mathrm{O}] comes from. Then:

• Since two hydrogen atoms get formed in this reaction, we get a factor of +2 in the first equation.

• Since one oxygen atom gets formed, we get a factor of +1 in the second equation.

• Since one water molecule gets used up, we get a factor of +1 in the third equation.

The dissociation of water (part 2)

Of course, we can also look at the reverse of the more realistic reaction involving a hydroxyl radical as an intermediate. Again, we just turn around the arrows in the Petri net we had:

Now the rate equation looks like this:

\begin{array}{ccl}  \frac{d [\mathrm{H}]}{d t} &=& + \alpha [\mathrm{OH}] + \beta [\mathrm{H}_2\mathrm{O}]  \\  \\  \frac{d [\mathrm{OH}]}{d t} &=& - \alpha [\mathrm{OH}] + \beta [\mathrm{H}_2 \mathrm{O}]  \\  \\  \frac{d [\mathrm{O}]}{d t} &=& + \alpha [\mathrm{OH}]  \\  \\  \frac{d [\mathrm{H}_2\mathrm{O}]}{d t} &=& - \beta [\mathrm{H}_2\mathrm{O}]  \end{array}

Do you see why? Test your understanding of the general recipe.

By the way: if you’re a category theorist, when I said “turn around all the arrows” you probably thought “opposite category”. And you’d be right! A Petri net is just a way of presenting a strict symmetric monoidal category that’s freely generated by some objects (the states) and some morphisms (the transitions). When we turn around all the arrows in our Petri net, we’re getting a presentation of the opposite symmetric monoidal category. For more details, try:

Vladimiro Sassone, On the category of Petri net computations, 6th International Conference on Theory and Practice of Software Development, Proceedings of TAPSOFT ’95, Lecture Notes in Computer Science 915, Springer, Berlin, pp. 334-348.

After I explain how stochastic Petri nets are related to quantum field theory, I hope to say more about this category theory business. But if you don’t understand it, don’t worry about it now—let’s do a few more examples.

The SI model

The SI model is an extremely simple model of an infectious disease. We can describe it using this Petri net:

There are two states: susceptible and infected. And there’s a transition called infection, where an infected person meets a susceptible person and infects them.

Suppose S is the number of susceptible people and I the number of infected ones. If the rate constant for infection is \beta, the rate equation is

\begin{array}{ccl}  \frac{d S}{d t} &=& - \beta S I \\ \\  \frac{d I}{d t} &=&  \beta S I   \end{array}

Do you see why?

By the way, it’s easy to solve these equations exactly. The total number of people doesn’t change, so S + I is a conserved quantity. Use this to get rid of one of the variables. You’ll get a version of the famous logistic equation, so the fraction of people infected must grow sort of like this:

Puzzle. Is there a stochastic Petri net with just one state whose rate equation is the logistic equation:

\frac{d P}{d t} = \alpha P - \beta P^2 \; ?

The SIR model

The SI model is just a warmup for the more interesting SIR model, which was invented by Kermack and McKendrick in 1927:

• W. O. Kermack and A. G. McKendrick, A Contribution to the mathematical theory of epidemics, Proc. Roy. Soc. Lond. A 115 (1927), 700-721.

The SIR model has an extra state, called resistant, and an extra transition, called recovery, where an infected person gets better and develops resistance to the disease:

If the rate constant for infection is \beta and the rate constant for recovery is \alpha, the rate equation for this stochastic Petri net is:

\begin{array}{ccl}  \frac{d S}{d t} &=& - \beta S I \\  \\  \frac{d I}{d t} &=&  \beta S I - \alpha I \\   \\  \frac{d R}{d t} &=&  \alpha I  \end{array}

See why?

I don’t know a ‘closed-form’ solution to these equations. But Kermack and McKendrick found an approximate solution in their original paper. They used this to model the death rate from bubonic plague during an outbreak in Bombay that started in 1905, and they got pretty good agreement. Nowadays, of course, we can solve these equations numerically on the computer.

The SIRS model

There’s an even more interesting model of infectious disease called the SIRS model. This has one more transition, called losing resistance, where a resistant person can go back to being susceptible. Here’s the Petri net:

Puzzle. If the rate constants for recovery, infection and loss of resistance are \alpha, \beta, and \gamma, write down the rate equations for this stochastic Petri net.

In the SIRS model we see something new: cyclic behavior! Say you start with a few infected people and a lot of susceptible ones. Then lots of people get infected… then lots get resistant… and then, much later, if you set the rate constants right, they lose their resistance and they’re ready to get sick all over again! You can sort of see it from the Petri net, which looks like a cycle.

I learned about the SI, SIR and SIRS models from this great book:

• Marc Mangel, The Theoretical Biologist’s Toolbox: Quantitative Methods for Ecology and Evolutionary Biology, Cambridge U. Press, Cambridge, 2006.

For more models of this type, see:

Compartmental models in epidemiology, Wikipedia.

A ‘compartmental model’ is closely related to a stochastic Petri net, but beware: the pictures in this article are not really Petri nets!

The general recipe revisited

Now let me remind you of the general recipe and polish it up a bit. So, suppose we have a stochastic Petri net with k states. Let x_i be the number of things in the ith state. Then the rate equation looks like:

\frac{d x_i}{d t} = ???

It’s really a bunch of equations, one for each 1 \le i \le k. But what is the right-hand side?

The right-hand side is a sum of terms, one for each transition in our Petri net. So, let’s assume our Petri net has just one transition! (If there are more, consider one at a time, and add up the results.)

Suppose the ith state appears as input to this transition m_i times, and as output n_i times. Then the rate equation is

\frac{d x_i}{d t} = r (n_i - m_i) x_1^{m_1} \cdots x_k^{m_k}

where r is the rate constant for this transition.

That’s really all there is to it! But subscripts make my eyes hurt more and more as I get older—this is the real reason for using index-free notation, despite any sophisticated rationales you may have heard—so let’s define a vector

x = (x_1, \dots , x_k)

that keeps track of how many things there are in each state. Similarly let’s make up an input vector:

m = (m_1, \dots, m_k)

and an output vector:

n = (n_1, \dots, n_k)

for our transition. And a bit more unconventionally, let’s define

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

Then we can write the rate equation for a single transition as

\frac{d x}{d t} = r (n-m) x^m

This looks a lot nicer!

Indeed, this emboldens me to consider a general stochastic Petri net with lots of transitions, each with their own rate constant. Let’s write T for the set of transitions and r(\tau) for the rate constant of the transition \tau \in T. Let n(\tau) and m(\tau) be the input and output vectors of the transition \tau. Then the rate equation for our stochastic Petri net is

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

That’s the fully general recipe in a nutshell. I’m not sure yet how helpful this notation will be, but it’s here whenever we want it.

Next time we’ll get to the really interesting part, where ideas from quantum theory enter the game! We’ll see how things in different states randomly transform into each other via the transitions in our Petri net. And someday we’ll check that the expected number of things in each state evolves according to the rate equation we just wrote down… at least in there limit where there are lots of things in each state.


Lifeboat Foundation

1 April, 2011

I’ve been invited to join this organization. But you can join too:

Lifeboat Foundation.

I hadn’t heard of it before. Do you know anything about it? Here’s their mission statement:

The Lifeboat Foundation is a nonprofit nongovernmental organization dedicated to encouraging scientific advancements while helping humanity survive existential risks and possible misuse of increasingly powerful technologies, including genetic engineering, nanotechnology, and robotics/AI, as we move towards the Singularity.

Lifeboat Foundation is pursuing a variety of options, including helping to accelerate the development of technologies to defend humanity, including new methods to combat viruses (such as RNA interference and new vaccine methods), effective nanotechnological defensive strategies, and even self-sustaining spacecolonies in case the other defensive strategies fail.

We believe that, in some situations, it might be feasible to relinquish technological capacity in the public interest (for example, we are against the U.S. government posting the recipe for the 1918 flu virus on the Internet).

We have some of the best minds on the planet working on programs to enable our survival. We invite you to join our cause!

It seems to have Nick Bostrom and Ray Kurzweil as two of its guiding figures: the overview features quotes from both.

Overview

An existential risk is a risk that is both global and terminal. Nick Bostrom defines it as a risk “where an adverse outcome would either annihilate Earth-originating intelligent life or permanently and drastically curtail its potential”. The term is frequently used to describe disaster and doomsday scenarios caused by non-friendly superintelligence, misuse of molecular nanotechnology, or other sources of danger.

The Lifeboat Foundation was formed to prevent existential events from happening, as once they occur, humanity may have no possibility to correct the error. Unfortunately governments, and humanity in general, always react AFTER a disaster has happened, and some disasters will leave no survivors so we must react BEFORE they occur. We must be proactive.

The Lifeboat Foundation is developing programs to prevent existential events (“shields”) as well as programs to preserve civilization (“preservers”) to survive such events.

Quotes

“Our approach to existential risks cannot be one of trial-and-error. There is no opportunity to learn from errors. The reactive approach — see what happens, limit damages, and learn from experience — is unworkable. Rather, we must take a proactive approach. This requires foresight to anticipate new types of threats and a willingness to take decisive preventive action and to bear the costs (moral and economic) of such actions.” — Nick Bostrom

“We cannot rely on trial-and-error approaches to deal with existential risks… We need to vastly increase our investment in developing specific defensive technologies… We are at the critical stage today for biotechnology, and we will reach the stage where we need to directly implement defensive technologies for nanotechnology during the late teen years of this century… A self-replicating pathogen, whether biological or nanotechnology based, could destroy our civilization in a matter of days or weeks.” — Ray Kurzweil

You’ll note there’s no mention here of global warming, mass extinction of species, oil depletion and other minor nuisances. Some people consider these problems insufficiently severe to count as “existential threats”… and thus, perhaps, best left to others. Some argue that there are already enough people worrying about these problem—while other threats need more attention than they’re getting.

That would be an interesting discussion to have. But I’m afraid there’s a cultural divide between the “green crowd” and the “tech crowd” that hinders such a discussion. The green crowd worries about things like global warming, the mass extinction that may currently be underway, and peak oil. The tech crowd worries about things like nanotechnology, artificial intelligence and asteroids hitting the Earth. Each crowd tends to think the other is a bit silly… and they don’t talk to each other enough. Am I just imagining this? I don’t think so.

Of course, any generalization this vast admits many exceptions. I like Gregory Benford because he confounds naive expectations: he thinks global warming is a desperately urgent problem that overshadows all others, but he’s willing to contemplate high-tech solutions. According to my theory, that should annoy both the green crowd and the tech crowd.

Personally I think all significant threats to civilization and biosphere should be evaluated and addressed in a unified way. Setting some aside because they’re “non-existential” or overly studied seems just as dangerous as setting others aside because they seem improbable or science-fiction-esque.

For one thing, I can imagine scenarios where medium-sized problems snowball into big “existential” ones. What’s the chance that in this century, global warming leads to droughts and famines which combined with oil shortages lead to political instability, the collapse of democratic governments, wars… and finally a world-wide nuclear or biological war? Maybe low… but I bet it’s higher than the chance of an asteroid hitting the Earth in this century.

I’m pleased to see that the Lifeboat Foundation plans “future programs” that will appeal to the green crowd:

ClimateShield
To protect against global warming and other unwanted climate changes.

BioPreserver
To preserve animal life and diversity on the planet.

EnergyPreserver
If our civilization ran out of energy, it would grind to a halt, so Lifeboat Foundation is looking for solutions.

However, their current programs are strongly focused on issues that appeal to the tech crowd. Maybe that’s okay, but maybe it’s a bit unbalanced:

AIShield
To protect against unfriendly AI (Artificial Intelligence).

AsteroidShield
To protect against devastating asteroid strikes.

BioShield
To protect against bioweapons and pandemics.

InternetShield
As the Internet grows in importance, an attack on it could cause physical as well as informational damage. An attack today on hospital systems or electric utilities could lead to deaths. In the future an attack could be used to alter the output that is produced by
nanofactories worldwide leading to massive deaths.

LifeShield Bunkers
Developing fallback positions on Earth in case programs such as our BioShield and NanoShield fail globally or locally.

NanoShield
To protect against ecophages and nonreplicating
nanoweapons.

ScientificFreedomShield
This shield strives to protect scientists from obstacles that would prevent latter day Max Plancks from completing their research.

SecurityPreserver
To prevent nuclear, biological, and nanotechnological attacks from occurring by using surveillance and
sousveillance
to identify terrorists before they are able to launch their attacks.

Space Habitats
To build fail-safes against global existential risks by encouraging the spread of sustainable human civilization beyond Earth.


Network Theory (Part 2)

31 March, 2011

Today I’d like to start telling you about some research Jacob Biamonte and I are doing on stochastic Petri nets, quantum field theory, and category theory. It’ll take a few blog entries to cover this story, which is part of a larger story about network theory.

Stochastic Petri nets are one of many different diagrammatic languages people have evolved to study complex systems. We’ll see how they’re used in chemistry, molecular biology, population biology and queuing theory, which is roughly the science of waiting in line. If you’re a faithful reader of this blog, you’ve already seen an example of a Petri net taken from chemistry:

It shows some chemicals and some reactions involving these chemicals. To make it into a stochastic Petri net, we’d just label each reaction by a nonnegative real number: the reaction rate constant, or rate constant for short.

In general, a Petri net will have a set of states, which we’ll draw as yellow circles, and a set of transitions, which we’ll draw as blue rectangles. Here’s a Petri net from population biology:

Now, instead of different chemicals, the states are different species. And instead of chemical reactions, the transitions are processes involving our species.

This Petri net has two states: rabbit and wolf. It has three transitions:

• In birth, one rabbit comes in and two go out. This is a caricature of reality: these bunnies reproduce asexually, splitting in two like amoebas.

• In predation, one wolf and one rabbit come in and two wolves go out. This is a caricature of how predators need to eat prey to reproduce.

• In death, one wolf comes in and nothing goes out. Note that we’re pretending rabbits don’t die unless they’re eaten by wolves.

If we labelled each transition with a rate constant, we’d have a stochastic Petri net.

To make this Petri net more realistic, we’d have to make it more complicated. I’m trying to explain general ideas here, not realistic models of specific situations. Nonetheless, this Petri net already leads to an interesting model of population dynamics: a special case of the so-called ‘Lotka-Volterra predator-prey model’. We’ll see the details soon.

More to the point, this Petri net illustrates some possibilities that our previous example neglected. Every transition has some ‘input’ states and some ‘output’ states. But a state can show up more than once as the output (or input) of some transition. And as we see in ‘death’, we can have a transition with no outputs (or inputs) at all.

But let me stop beating around the bush, and give you the formal definitions. They’re simple enough:

Definition. A Petri net consists of a set S of states and a set T of transitions, together with a function

i : S \times T \to \mathbb{N}

saying how many copies of each state shows up as input for each transition, and a function

o: S \times T \to \mathbb{N}

saying how many times it shows up as output.

Definition. A stochastic Petri net is a Petri net together with a function

r: T \to [0,\infty)

giving a rate constant for each transition.

Starting from any stochastic Petri net, we can get two things. First:

• The master equation. This says how the probability that we have a given number of things in each state changes with time.

Since stochastic means ‘random’, the master equation is what gives stochastic Petri nets their name. The master equation is the main thing I’ll be talking about in future blog entries. But not right away!

Why not?

In chemistry, we typically have a huge number of things in each state. For example, a gram of water contains about 3 \times 10^{22} water molecules, and a smaller but still enormous number of hydroxide ions (OH), hydronium ions (H3O+), and other scarier things. These things blunder around randomly, bump into each other, and sometimes react and turn into other things. There’s a stochastic Petri net describing all this, as we’ll eventually see. But in this situation, we don’t usually want to know the probability that there are, say, exactly 310,1849,578,476,264 hydronium ions. That would be too much information! We’d be quite happy knowing the expected value of the number of hydronium ions, so we’d be delighted to have a differential equation that says how this changes with time.

And luckily, such an equation exists—and it’s much simpler than the master equation. So, today we’ll talk about:

• The rate equation. This says how the expected number of things in each state changes with time.

But first, I hope you get the overall idea. The master equation is stochastic: at each time the number of things in each state is a random variable taking values in \mathbb{N}, the set of natural numbers. The rate equation is deterministic: at each time the expected number of things in each state is a non-random variable taking values in [0,\infty), the set of nonnegative real numbers. If the master equation is the true story, the rate equation is only approximately true—but the approximation becomes good in some limit where the expected value of the number of things in each state is large, and the standard deviation is comparatively small.

If you’ve studied physics, this should remind you of other things. The master equation should remind you of the quantum harmonic oscillator, where energy levels are discrete, and probabilities are involved. The rate equation should remind you of the classical harmonic oscillator, where energy levels are continuous, and everything is deterministic.

When we get to the ‘original research’ part of our story, we’ll see this analogy is fairly precise! We’ll take a bunch of ideas from quantum mechanics and quantum field theory, and tweak them a bit, and show how we can use them to describe the master equation for a stochastic Petri net.

Indeed, the random processes that the master equation describes can be drawn as pictures:

This looks like a Feynman diagram, with animals instead of particles! It’s pretty funny, but the resemblance is no joke: the math will back it up.

I’m dying to explain all the details. But just as classical field theory is easier than quantum field theory, the rate equation is simpler than the master equation. So we should start there.

The rate equation

If you hand me a stochastic Petri net, I can write down its rate equation. Instead of telling you the general rule, which sounds rather complicated at first, let me do an example. Take the Petri net we were just looking at:

We can make it into a stochastic Petri net by choosing a number for each transition:

• the birth rate constant \beta

• the predation rate constant \gamma

• the death rate constant \delta

Let x(t) be the number of rabbits and let y(t) be the number of wolves at time t. Then the rate equation looks like this:

\frac{d x}{d t} = \beta x - \gamma x y

\frac{d y}{d t} = \gamma x y - \delta y

It’s really a system of equations, but I’ll call the whole thing “the rate equation” because later we may get smart and write it as a single equation.

See how it works?

• We get a term \beta x in the equation for rabbits, because rabbits are born at a rate equal to the number of rabbits times the birth rate constant \beta.

• We get a term - \delta y in the equation for wolves, because wolves die at a rate equal to the number of wolves times the death rate constant \delta.

• We get a term - \gamma x y in the equation for rabbits, because rabbits die at a rate equal to the number of rabbits times the number of wolves times the predation rate constant \gamma.

• We also get a term \gamma x y in the equation for wolves, because wolves are born at a rate equal to the number of rabbits times the number of wolves times \gamma.

Of course I’m not claiming that this rate equation makes any sense biologically! For example, think about predation. The \gamma x y terms in the above equation would make sense if rabbits and wolves roamed around randomly, and whenever a wolf and a rabbit came within a certain distance, the wolf had a certain probability of eating the rabbit and giving birth to another wolf. At least it would be make sense in the limit of large numbers of rabbits and wolves, where we can treat x and y as varying continuously rather than discretely. That’s a reasonable approximation to make sometimes. Unfortunately, rabbits and wolves don’t roam around randomly, and a wolf doesn’t spit out a new wolf each time it eats a rabbit.

Despite that, the equations

\frac{d x}{d t} = \beta x - \gamma x y

\frac{d y}{d t} = \gamma x y - \delta y

are actually studied in population biology. As I said, they’re a special case of the Lotka-Volterra predator-prey model, which looks like this:

\frac{d x}{d t} = \beta x - \gamma x y

\frac{d y}{d t} = \epsilon x y - \delta y

The point is that while these models are hideously oversimplified and thus quantitatively inaccurate, they exhibit interesting qualititative behavior that’s fairly robust. Depending on the rate constants, these equations can show either a stable equilibrium or stable periodic behavior. And we go from one regime to another, we see a kind of catastrophe called a “Hopf bifurcation”. I explained all this in week308 and week309. There I was looking at some other equations, not the Lotka-Volterra equations. But their qualitative behavior is the same!

If you want stochastic Petri nets that give quantitatively accurate models, it’s better to retreat to chemistry. Compared to animals, molecules come a lot closer to roaming around randomly and having a chance of reacting when they come within a certain distance. So in chemistry, rate equations can be used to make accurate predictions.

But I’m digressing. I should be explaining the general recipe for getting a rate equation from a stochastic Petri net! You might not be able to guess it from just one example. But I sense that you’re getting tired. So let’s stop now. Next time I’ll do more examples, and maybe even write down a general formula. But if you’re feeling ambitious, you can try this now:

Puzzle. Can you write down a stochastic Petri net whose rate equation is the Lotka-Volterra predator-prey model:

\frac{d x}{d t} = \beta x - \gamma x y

\frac{d y}{d t} = \epsilon x y - \delta y

for arbitrary \beta, \gamma, \delta, \epsilon \ge 0? If not, for which values of these rate constants can you do it?

References

If you want to study a bit on your own, here are some great online references on stochastic Petri nets and their rate equations:

• Peter J. E. Goss and Jean Peccoud, Quantitative modeling of stochastic systems in molecular biology by using stochastic Petri nets, Proc. Natl. Acad. Sci. USA 95 (June 1998), 6750-6755.

• Jeremy Gunawardena, Chemical reaction network theory for in-silico biologists.

• Martin Feinberg, Lectures on reaction networks.

I should admit that the first two talk about ‘chemical reaction networks’ instead of Petri nets. That’s no big deal: any chemical reaction network gives a Petri net in a pretty obvious way. You can probably figure out how; if you get stuck, just ask.

Also, I should admit that Petri net people say place where I’m saying state.

Here are some other references, which aren’t free unless you have an online subscription or access to a library:

• Peter J. Haas, Stochastic Petri Nets: Modelling, Stability, Simulation, Springer, Berlin, 2002.

• F. Horn and R. Jackson, General mass action kinetics, Archive for Rational Mechanics and Analysis 47 (1972), 81–116.

• Ina Koch, Petri nets – a mathematical formalism to analyze chemical reaction networks, Molecular Informatics 29 (2010), 838–843.

• Darren James Wilkinson, Stochastic Modelling for Systems Biology, Taylor & Francis, New York, 2006.


Environment and Sustainability Institute

30 March, 2011

My friend John Barrett pointed out that the Environmental and Sustainability Institute at Exeter has jobs for mathematicians and statisticians who “combine research expertise in areas such as computational statistics, data modelling, system dynamics, control, optimization and/or computation with vision and innovation so as to transform research across the environment and sustainability agenda”.

I got his email today; the deadline has already passed, but maybe you can slip under the wire:

Professor in Applied Mathematics/Statistics

Ref: R10168

Deadline: 21/03/2011

Location: Environment and Sustainability Institute, College of Engineering, Mathematics and Physical Sciences, University of Exeter, Tremough Campus (Cornwall).

Salary: by negotiation, generous holiday allowances, flexible working, pension scheme, car lease scheme and relocation package.

We are taking a unique opportunity to create a world-leading, interdisciplinary centre investigating the consequences of environmental change and the mitigation and management of its effects. Building on the existing academic strengths of the University of Exeter, the 30 million pound Environment and Sustainability Institute (ESI) will focus on undertaking significant cutting-edge research in finding solutions to problems of environmental change, and the application of an ecosystem services approach. There will be three major themes: Clean technologies, Natural environment, and Social science. The initiative is funded with a 22.9 million pound investment from the European Regional Development Fund and 6.6 million pounds from the South West of England Regional Development Agency.

We are looking for research leaders in Applied Mathematics and/or Statistics who will combine research expertise in areas such as computational statistics, data modelling, system dynamics, control, optimization and/or computation with vision and innovation so as to transform research across the environment and sustainability agenda. Applicants will have a strong track record of research funding and international quality publications, together with proven ability in teaching and curriculum development.

Applicants are encouraged to contact the Dean of the College, Prof Ken Evans to discuss the posts further. Informal enquiries can be made to Prof Peter Ashwin You may also wish to consult our web site for further details of the College.

The full range of necessary skills and experience can be found in the Job Description and Person Specification document.

Your full academic CV should be accompanied by a short letter of application explaining why you are interested in the post.

For further information about the Environment and Sustainability Institute visit our website.

The University of Exeter is an equal opportunity employer which is ‘Positive About Disabled People’. Whilst all applicants will be judged on merit alone, we particularly welcome applications from groups currently under-represented in the workforce.


Follow

Get every new post delivered to your Inbox.

Join 3,095 other followers