Dysnomia

24 March, 2012

This morning I was looking at a nice website that lets you zoom in or out and see objects of different sizes, ranging from the Planck length to the entire observable Universe:

• Cary Huang, The Scale of the Universe.

When I zoomed out to a size somewhat larger than Rhode Island, I saw something strange:

Huh? I recently saw a movie called Melancholia, about a pair of sisters struggling with fear and depression, and a blue planet of that name that’s headed for Earth. And indeed, dysnomia is another mental disorder: a difficulty in remembering words, names or numbers. Are astronomical bodies named for mental disorders catching on?

I don’t know—but unlike Melancholia, Dysnomia is real! It’s the speck at left here:

The larger blob of light is the dwarf planet Eris. Dysnomia is its moon. Both were discovered in 2005 by a team at Palomar led by Mike Brown.

Why the funny name? In Greek mythology, Dysnomia (Δυσνομία), meaning ‘lawlessness’, was the daughter of Eris, the goddess of strife and discord. You may remember the dwarf planet Eris under its tentative early name ‘Xena’, from the TV character. That was deemed too silly to stand.

Eris is 30% more massive than Pluto, and thus it helped lead to a redefinition of ‘planet’: both are now called dwarf planets, because they aren’t big enough to clear their neighborhood of debris.

Eris has a highly eccentric orbit, and it takes 580 years to go around the Sun. Right now it’s at its farthest from the Sun: three times as far as Pluto. So it’s very cold, about 30 kelvin (-243 °C), and its surface is covered with methane ice. But in about 290 years its temperature will nearly double, soaring to a balmy 56 kelvin (-217 °C). Its methane ice will melt and then evaporate away, giving this world a new atmosphere! That’s pretty amazing: an annual atmosphere.

The surface of Eris is light gray. So, it’s quite different than Pluto and Neptune’s moon Triton, which are reddish due to deposits of tholins—complex compounds formed by bombarding hydrocarbons and nitrogen with sunlight.

Remember the smoggy orange surface of Titan?

That’s due to tholins! It’s possible that on Eris the tholins are currently covered by methane ice.

I wish I knew more about what tholins actually are—their actual chemical structures. But they’re a complicated mess of stuff, and they vary from place to place: people talk about Titan tholins, Triton tholins and ice tholins.

Indeed, the term “tholin” comes from the Greek word tholós (θολός), meaning “not clear”. It was coined by Carl Sagan to describe the mysterious mix of substances he created in experiments on the gas mixtures found in Titan’s atmosphere… experiments a bit like the old Urey-Miller experiment which created amino acids from stuff that was present on the early Earth—water, methane, ammonia, hydrogen—together with lots of electrical discharges.

Might tholins become complicated enough to lead to life in the outer Solar System, at least on relatively peppy worlds like Titan? There’s a complex cycle going on there:

Here ‘Da’ means daltons: a dalton is a unit of mass equal to about one hydrogen atom, so a molecule that’s 2000 Da could be made of 2000 hydrogens or, more reasonably, 2000/12 ≈ 166 carbons, or various other things. The point is: that’s a pretty big complicated molecule—big enough to be very interesting!

On Earth, many soil bacteria are able to use tholins as their sole source of carbon. Some think that tholins may have been the first food for microbes! In fact, some scientists have speculated that Earth may have been seeded by tholins from comets early in its development, providing the raw material necessary for life. But ever since the Great Oxygenation Event, the Earth’s surface has been too oxidizing for tholins to exist here.

Tholins have also been found outside our Solar System:

Red dust in disk may harbor precursors to life, Carnegie Institute news release, 5 January 2008.

I bet tholins often go hand-in-hand with PAHs, or polycyclic aromatic hydrocarbons. PAHs are also common in outer space. In Earth you can find them in soot, or the tarry stuff that forms in a barbecue grill. Wherever carbon-containing materials suffer incomplete combustion, you’ll find PAHs.

PAHs are made of hexagonal rings of carbon atoms, with some hydrogens along the edges:

Benzene has a single hexagonal ring, with 6 carbons and 6 hydrogens. You’ve probably heard about naphthalene, which is used for mothballs: this consists of two hexagonal rings stuck together. True PAHs have more. With three rings you can make anthracene:

and phenanthrene:

With four, you can make napthacene:

pyrene:

triphenylene:

and chrysene:

And so on! The game just gets more complicated as you get to use more puzzle pieces.

In 2004, a team of scientists discovered anthracene and pyrene in an amazing structure called the Red Rectangle!

Here two stars 2300 light years from us are spinning around each other while pumping out a huge torus of icy dust grains and hydrocarbon molecules. It’s not really shaped like a rectangle or X—it just looks that way from here. The whole scene is about 1/3 of a light year across.

This was first time such complex molecules had been
found in space:

• Uma P. Vijh, Adolf N. Witt, and Karl D. Gordon, Small polycyclic aromatic hydrocarbons in the Red Rectangle, The Astrophysical Journal 619 (2005), 368-378.

By now, lots of organic molecules have been found in interstellar or circumstellar space. There’s a whole "ecology" of organic chemicals out there, engaged in complex reactions. Life on planets might someday be seen as just an aspect of this larger ecology. PAHs and tholins are probably among the dominant players in this ecology, at least at this stage.

Indeed, I’ve read that about 10% of the interstellar carbon is in the form of PAHs—big ones, with about 50 carbons per molecule. They’re common because they’re incredibly stable. They’ve even been found riding the shock wave of a supernova explosion!

PAHs are also found in meteorites called "carbonaceous chondrites". These space rocks contain just a little carbon: about 3% by weight. But 80% of this carbon is in the form of PAHs.

Here’s an interview with a scientist who thinks PAHs were important precursors of life on Earth:

Aromatic world, interview with Pascale Ehrenfreund, Astrobiology Magazine.

Also try this:

PAH world hypothesis, Wikipedia.

This speculative hypothesis says that PAHs were abundant in the primordial soup of the early Earth and played a major role in the origin of life by mediating the synthesis of RNA molecules, leading to the (also speculative) RNA world.

Another radical theory has been proposed by Prof. Sun Kwok, author of Organic Matter in the Universe. He claims that instead of PAHs, complex molecules like this would do better at explaining the spectra of interstellar clouds:

Would this molecule count as a tholin? Maybe so: I don’t know. He says:

Our work has shown that stars have no problem making complex organic compounds under near-vacuum conditions. Theoretically, this is impossible, but observationally we can see it happening.

For more see:

• Sun Kwok and Yong Zhang, Mixed aromatic–aliphatic organic nanoparticles as carriers of unidentified infrared emission features, Nature 479 (2011), 80–83.

This paper isn’t free, but here’s a summary that is:

Astronomers discover complex organic matter exists throughout the Universe, ScienceDaily, 26 October 2011.

However, I’d take this with a grain of salt until more confirmation comes along! They’re matching very complicated spectra to hypothetical chemicals, without yet any understanding of how these chemicals could be formed in space. It would be very cool if true.

Regardless of how the details play out, I think we’ll eventually see that organic life across the universe is a natural outgrowth of the organic chemistry of PAHs, tholins and related chemicals. It will be great to see the whole story: how much in common life has in different locations, and how much variation there is. It may be rare, but the universe is very large, so there must be statistical patterns in how life works.

It goes to show how everything is connected. Starting from a chance encounter with Dysnomia, we’ve been led to ponder another planet whose atmosphere liquifies and then freezes every year… and then wonder about why so many objects in the outer solar system are red… and why the same chemicals you find in the tarry buildup on a barbecue grill are also seen in outer space… and whether life on Earth could have been started by complex compounds from comets… and whether life on planets is just part of a larger interstellar chemical ‘ecology’. Not bad for a Saturday morning!


Network Theory (Part 18)

21 November, 2011

joint with Jacob Biamonte

Last time we explained how ‘reaction networks’, as used in chemistry, are just another way of talking about Petri nets. We stated an amazing result, the deficiency zero theorem, which settles quite a number of questions about chemical reactions. Now let’s illustrate it with an example.

Our example won’t show how powerful this theorem is: it’s too simple. But it’ll help explain the ideas involved.

Diatomic molecules

A diatomic molecule consists of two atoms of the same kind, stuck together:

At room temperature there are 5 elements that are diatomic gases: hydrogen, nitrogen, oxygen, fluorine, chlorine. Bromine is a diatomic liquid, but easily evaporates into a diatomic gas:

Iodine is a crystal at room temperatures:

but if you heat it a bit, it becomes a diatomic liquid and then a gas:

so people often list it as a seventh member of the diatomic club.

When you heat any diatomic gas enough, it starts becoming a ‘monatomic’ gas as molecules break down into individual atoms. However, just as a diatomic molecule can break apart into two atoms:

A_2 \to A + A

two atoms can recombine to form a diatomic molecule:

A + A \to A_2

So in equilibrium, the gas will be a mixture of diatomic and monatomic forms. The exact amount of each will depend on the temperature and pressure, since these affect the likelihood that two colliding atoms stick together, or a diatomic molecule splits apart. The detailed nature of our gas also matters, of course.

But we don’t need to get into these details here! Instead, we can just write down the ‘rate equation’ for the reactions we’re talking about. All the details we’re ignoring will be hiding in some constants called ‘rate constants’. We won’t try to compute these; we’ll leave that to our chemist friends.

A reaction network

To write down our rate equation, we start by drawing a ‘reaction network’. For this, we can be a bit abstract and call the diatomic molecule B instead of A_2. Then it looks like this:

We could write down the same information using a Petri net:

But today let’s focus on the reaction network! Staring at this picture, we can read off various things:

Species. The species are the different kinds of atoms, molecules, etc. In our example the set of species is S = \{A, B\}.

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

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

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

So, in our example we have:

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

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

• Reactions: T =  \{A+A\to B, B\to A+A\}.

To get the rate equation, we also need one more piece of information: a rate constant r(\tau) for each reaction \tau \in T. This is a nonnegative real number that affects how fast the reaction goes. All the details of how our particular diatomic gas behaves at a given temperature and pressure are packed into these constants!

The rate equation

The rate equation says how the expected numbers of the various species—atoms, molecules and the like—changes with time. This equation is deterministic. It’s a good approximation when the numbers are large and any fluctuations in these numbers are negligible by comparison.

Here’s the general form of the rate equation:

\displaystyle{ \frac{d}{d t} x_i =  \sum_{\tau\in T} r(\tau) \, (n_i(\tau)-m_i(\tau)) \, x^{m(\tau)}  }

Let’s take a closer look. The quantity x_i is the expected population of the ith species. So, this equation tells us how that changes. But what about the right hand side? As you might expect, it’s a sum over reactions. And:

• The term for the reaction \tau is proportional to the rate constant r(\tau).

• Each reaction \tau goes between two complexes, so we can write it as m(\tau) \to n(\tau). Among chemists the input m(\tau) is called the reactant complex, and the output is called the product complex. The difference n_i(\tau)-m_i(\tau) tells us how many items of species i get created, minus how many get destroyed. So, it’s the net amount of this species that gets produced by the reaction \tau. The term for the reaction \tau is proportional to this, too.

• Finally, the law of mass action says that the rate of a reaction is proportional to the product of the concentrations of the species that enter as inputs. More precisely, if we have a reaction \tau with input is the complex m(\tau), we define x^{m(\tau)} = x_1^{m_1(\tau)} \cdots x_k^{m_k(\tau)}. The law of mass action says the term for the reaction \tau is proportional to this, too!

Let’s see what this says for the reaction network we’re studying:

Let’s write x_1(t) for the number of A atoms and x_2(t) for the number of B molecules. Let the rate constant for the reaction B \to A + A be \alpha, and let the rate constant for A + A \to B be \beta, Then the rate equation is this:

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

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

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

We’ve got two species, A and B. But remember, B is just an abbreviation for a molecule made of two A atoms. So, the total number of A atoms is conserved by the reactions in our network. This is the number of A‘s plus twice the number of B‘s: x_1 + 2x_2. So, this should be a conserved quantity: it should not change with time. Indeed, by adding the first equation above to twice the second, we see:

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

As a consequence, any solution will stay on a line

x_1 + 2 x_2 = c

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

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

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

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

and then solve for x_1.

This sort of trick won’t work for more complicated examples.
But the idea remains important: the numbers of atoms of various kinds—hydrogen, helium, lithium, and so on—are conserved by chemical reactions, so a solution of the rate equation can’t roam freely in \mathbb{R}^S. It will be trapped on some hypersurface, which is called a ‘stoichiometric compatibility class’. And this is very important.

We don’t feel like doing the integral required to solve our rate equation in closed form, because this idea doesn’t generalize too much. On the other hand, we can always solve the rate equation numerically. So let’s try that!

For example, suppose we set \alpha = \beta = 1. We can plot the solutions for three different choices of initial conditions, say (x_1,x_2) = (0,3), (4,0), and (3,3). We get these graphs:

It looks like the solution always approaches an equilibrium. We seem to be getting different equilibria for different initial conditions, and the pattern is a bit mysterious. However, something nice happens when we plot the ratio x_1^2 / x_2:

Apparently it always converges to 1. Why should that be? It’s not terribly surprising. With both rate constants equal to 1, the reaction A + A \to B proceeds at a rate equal to the square of the number of A‘s, namely x_1^2. The reverse reaction proceeds at a rate equal to the number of B‘s, namely x_2. So in equilibrium, we should have x_1^2 = x_2.

But why is the equilibrium stable? In this example we could see that using the closed-form solution, or maybe just common sense. But it also follows from a powerful theorem that handles a lot of reaction networks.

The deficiency zero theorem

It’s called Feinberg’s deficiency zero theorem, and we saw it last time. Very roughly, it says that if our reaction network is ‘weakly reversible’ and has ‘deficiency zero’, the rate equation will have equilibrium solutions that behave about as nicely as you could want.

Let’s see how this works. We need to remember some jargon:

Weakly reversible. A reaction network is weakly reversible if for every reaction v \to w in the network, there exists a path of reactions in the network starting at w and leading back to v.

Reversible. A reaction network is reversible if for every reaction v \to w in the network, w \to v is also a reaction in the network. Any reversible reaction network is weakly reversible. Our example is reversible, since it consists of reactions A + A \to B, B \to A + A.

But what about ‘deficiency zero’? We defined that concept last time, but let’s review:

Connected component. A reaction network gives a kind of graph with complexes as vertices and reactions as edges. Two complexes lie in the same connected component if we can get from one to the other by a path of reactions, where at each step we’re allowed to go either forward or backward along a reaction. Chemists call a connected component a linkage class. In our example there’s just one:

Stoichiometric subspace. The stoichiometric subspace is the subspace \mathrm{Stoch} \subseteq \mathbb{R}^S spanned by the vectors of the form w - v for all reactions v \to w in our reaction network. This subspace describes the directions in which a solution of the rate equation can move. In our example, it’s spanned by B - 2 A and 2 A - B, or if you prefer, (-2,1) and (2,-1). These vectors are linearly dependent, so the stoichiometric subspace has dimension 1.

Deficiency. The deficiency of a reaction network is the number of complexes, minus the number of connected components, minus the dimension of the stoichiometric subspace. In our example there are 2 complexes, 1 connected component, and the dimension of the stoichiometric subspace is 1. So, our reaction network has deficiency 2 – 1 – 1 = 0.

So, the deficiency zero theorem applies! What does it say? To understand it, we need a bit more jargon. First of all, a vector x \in \mathbb{R}^S tells us how much we’ve got of each species: the amount of species i \in S is the number x_i. And then:

Stoichiometric compatibility class. Given a vector v\in \mathbb{R}^S, its stoichiometric compatibility class is the subset of all vectors that we could reach using the reactions in our reaction network:

\{ v + w \; : \; w \in \mathrm{Stoch} \}

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

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

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

x_1 + 2x_2 = c

We’ve already seen that if we start with initial conditions on such a line, the solution will stay on this line. And that’s how it always works: as time passes, any solution of the rate equation stays in the same stoichiometric compatibility class!

In other words: the stoichiometric subspace is defined by a bunch of linear equations, one for each linear conservation law that all the reactions in our network obey.

Here a linear conservation law is a law saying that some linear combination of the numbers of species does not change.

Next:

Positivity. A vector in \mathbb{R}^S is positive if all its components are positive; this describes a a container of chemicals where all the species are actually present. The positive stoichiometric compatibility class of x\in \mathbb{R}^S consists of all positive vectors in its stoichiometric compatibility class.

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

Zero Deficiency Theorem (Feinberg). If a reaction network has deficiency zero and is weakly reversible, and the rate constants are positive, the rate equation has exactly one equilibrium solution in each positive stoichiometric compatibility class. Any sufficiently nearby solution that starts in the same class will approach this equilibrium as t \to +\infty.

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

x_1 + 2x_2 = c

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

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

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

Solving these, we get

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

So, these are our equilibrium solutions. It’s easy to verify that indeed, there’s one of these in each stoichiometric compatibility class x_1 + 2x_2 = c. And the zero deficiency theorem also tells us that any sufficiently nearby solution that starts in the same class will approach this equilibrium as t \to \infty.

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

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

will actually have

\displaystyle{\lim_{t \to +\infty} \frac{x_1^2}{x_2} = 1 }

But in fact, in this example we don’t even need to start near the equilibrium for our solution to approach the equilibrium! What about in general? We don’t know, but just to get the ball rolling, we’ll risk the following wild guess:

Conjecture. If a reaction network has deficiency zero and is weakly reversible, and the rate constants are positive, the rate equation has exactly one equilibrium solution in each positive stoichiometric compatibility class, and any positive solution that starts in the same class will approach this equilibrium as t \to +\infty.

If anyone knows a proof or counterexample, we’d be interested. If this result were true, it would really clarify the dynamics of reaction networks in the zero deficiency case.

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


A Math Puzzle Coming From Chemistry

23 October, 2011

I posed this puzzle a while back, and nobody solved it. That’s okay—now that I think about it, I’m not sure how to solve it either!

It seems to involve group theory. But instead of working on it, solving it and telling you the answer, I’d rather dump all the clues in your lap, so we can figure it out together.

Suppose we have an ethyl cation. We’ll pretend it looks like this:

As I explained before, it actually doesn’t—not in real life. But never mind! Realism should never stand in the way of a good puzzle.

Continuing on in this unrealistic vein, we’ll pretend that the two black carbon atoms are distinguishable, and so are the five white hydrogen atoms. As you can see, 2 of the hydrogens are bonded to one carbon, and 3 to the other. We don’t care how the hydrogens are arranged, apart from which carbon each hydrogen is attached to. Given this, there are

2 \times \displaystyle{ \binom{5}{2} = 20 }

ways to arrange the hydrogens. Let’s call these arrangements states.

Now draw a dot for each of these 20 states. Draw an edge connecting two dots whenever you can get from one state to another by having a hydrogen hop from the carbon with 2 hydrogens to the carbon with 3. You’ll get this picture, called the Desargues graph:

The red dots are states where the first carbon has 2 hydrogens attached to it; the blue ones are states where the second carbon has 2 hydrogens attached to it. So, each edge goes between a red and a blue dot. And there are 3 edges coming out of each dot, since there are 3 hydrogens that can make the jump!

Now, the puzzle is to show that you can also get the Desargues graph from a different kind of molecule. Any molecule shaped like this will do:

The 2 balls on top and bottom are called axial, while the 3 around the middle are called equatorial.

There are various molecules like this. For example, phosphorus pentachloride. Let’s use that.

Like the ethyl cation, phosphorus pentachloride also has 20 states… but only if count them a certain way! We have to treat all 5 chlorines as distinguishable, but think of two arrangements of them as the same if we can rotate one to get the other. Again, I’m not claiming this is physically realistic: it’s just for the sake of the puzzle.

Phosphorus pentachloride has 6 rotational symmetries, since you can turn it around its axis 3 ways, but also flip it over. So, it has

\displaystyle{ \frac{5!}{6}  = 20}

states.

That’s good: exactly the number of dots in the Desargues graph! But how about the edges? We get these from certain transitions between states. These transitions are called pseudorotations, and they look like this:

Phosphorus pentachloride really does this! First the 2 axial guys move towards each other to become equatorial. Beware: now the equatorial ones are no longer in the horizontal plane: they’re in the plane facing us. Then 2 of the 3 equatorial guys swing out to become axial.

To get from one state to another this way, we have to pick 2 of the 3 equatorial guys to swing out and become axial. There are 3 choices here. So, we again get a graph with 20 vertices and 3 edges coming out of each vertex.

Puzzle. Is this graph the Desargues graph? If so, show it is.

I read in some chemistry papers that it is. But is it really? And if so, why? David Corfield suggested a promising strategy. He pointed out that we just need to get a 1-1 correspondence between

states of the ethyl cation and states of phosphorus pentachloride,

together with a compatible 1-1 correspondence between

transitions of the ethyl cation and transitions of phosphorus pentachloride.

And he suggested that to do this, we should think of the split of hydrogens into a bunch of 2 and a bunch of 3 as analogous to the split of chlorines into a bunch of 2 (the ‘axial’ ones) and a bunch of 3 (the ‘equatorial’ ones).

It’s a promising idea. There’s a problem, though! In the ethyl cation, a single hydrogen hops from the bunch of 3 to the bunch of 2. But in a pseudorotation, two chlorines go from the bunch of 2 to the bunch of 3… and meanwhile, two go back from the bunch of 3 to bunch of 2.

And if you think about it, there’s another problem too. In the ethyl cation, there are 2 distinguishable carbons. One of them has 3 hydrogens attached, and one doesn’t. But in phosphorus pentachloride it’s not like that. The 3 equatorial chlorines are just that: equatorial. They don’t have 2 choices about how to be that way. Or do they?

Well, there’s more to say, but this should already make it clear that getting ‘natural’ one-to-one correspondences is a bit tricky… if it’s even possible at all!

If you know some group theory, we could try solving the problem using the ideas behind Felix Klein’s ‘Erlangen program’. The group of permutations of 5 things, say S_5, acts as symmetries of either molecule. For the ethyl cation the set of states will be X  = S_5/G for some subgroup G. You can think of X as a set of structures of some sort on a 5-element set. The group S_5 acts on X, and the transitions will give an invariant binary relation on X, For phosphorus pentachloride we’ll have some set of states X' = S_5/G' for some other subgroup G', and the transitions will give an invariant relation on X'.

We could start by trying to see if G is the same as G'—or more precisely, conjugate. If they are, that’s a good sign. If not, it’s bad: it probably means there’s no ‘natural’ way to show the graph for phosphorus pentachloride is the Desargues graph.

I could say more, but I’ll stop here. In case you’re wondering, all this is just a trick to get more mathematicians interested in chemistry. A few may then go on to do useful things.


Network Theory (Part 14)

15 October, 2011

We’ve been doing a lot of hard work lately. Let’s take a break and think about a fun example from chemistry!

The ethyl cation

Suppose you start with a molecule of ethane, which has 2 carbons and 6 hydrogens arranged like this:

Then suppose you remove one hydrogen. The result is a positively charged ion, or ‘cation’. When I was a kid, I thought the opposite of a cation should be called a ‘dogion’. Alas, it’s not.

This particular cation, formed from removing one hydrogen from an ethane molecule, is called an ‘ethyl cation’. People used to think it looked like this:

They also thought a hydrogen could hop from the carbon with 3 hydrogens attached to it to the carbon with 2. So, they drew a graph with a vertex for each way the hydrogens could be arranged, and an edge for each hop. It looks really cool:

The red vertices come from arrangements where the first carbon has 2 hydrogens attached to it, and the blue vertices come from those where the second carbon has 2 hydrogens attached to it. So, each edge goes between a red vertex and a blue vertex.

This graph has 20 vertices, which are arrangements or ‘states’ of the ethyl cation. It has 30 edges, which are hops or ‘transitions’. Let’s see why those numbers are right.

First I need to explain the rules of the game. The rules say that the 2 carbon atoms are distinguishable: there’s a ‘first’ one and a ‘second’ one. The 5 hydrogen atoms are also distinguishable. But, all we care about is which carbon atom each hydrogen is bonded to: we don’t care about further details of its location. And we require that 2 of the hydrogens are bonded to one carbon, and 3 to the other.

If you’re a physicist, you may wonder why the rules work this way: after all, at a fundamental level, identical particles aren’t really distinguishable. I’m afraid I can’t give a fully convincing explanation right now: I’m just reporting the rules as they were told to me!

Given these rules, there are 2 choices of which carbon has two hydrogens attached to it. Then there are

\displaystyle{ \binom{5}{2} = \frac{5 \times 4}{2 \times 1} = 10}

choices of which two hydrogens are attached to it. This gives a total of 2 × 10 = 20 states. These are the vertices of our graph: 10 red and 10 blue.

The edges of the graph are transitions between states. Any hydrogen in the group of 3 can hop over to the group of 2. There are 3 choices for which hydrogen atom makes the jump. So, starting from any vertex in the graph there are 3 edges. This means there are 3 \times 20 / 2 = 30 edges.

Why divide by 2? Because each edge touches 2 vertices. We have to avoid double-counting them.

The Desargues graph

The idea of using this graph in chemistry goes back to this paper:

• A. T. Balaban, D. Fǎrcaşiu and R. Bǎnicǎ, Graphs of multiple 1,2-shifts in carbonium ions and related systems, Rev. Roum. Chim. 11 (1966), 1205.

This paper is famous because it was the first to use graphs in chemistry to describe molecular transitions, as opposed to using them as pictures of molecules!

But this particular graph was already famous for other reasons. It’s called the Desargues-Levi graph, or Desargues graph for short:

Desargues graph, Wikipedia.

Later I’ll say why it’s called this.

There are lots of nice ways to draw the Desargues graph. For example:

The reason why we can draw such pretty pictures is that the Desargues graph is very symmetrical. Clearly any permutation of the 5 hydrogens acts as a symmetry of the graph, and so does any permutation of the 2 carbons. This gives a symmetry group S_5 \times S_2, which has 5! \times 2! = 240 elements. And in fact this turns out to be the full symmetry group of the Desargues graph.

The Desargues graph, its symmetry group, and its applications to chemistry are discussed here:

• Milan Randic, Symmetry properties of graphs of interest in chemistry: II: Desargues-Levi graph, Int. Jour. Quantum Chem. 15 (1997), 663-682.

The ethyl cation, revisited

We can try to describe the ethyl cation using probability theory. If at any moment its state corresponds to some vertex of the Desargues graph, and it hops randomly along edges as time goes by, it will trace out a random walk on the Desargues graph. This is a nice example of a Markov process!

We could also try to describe the ethyl cation using quantum mechanics. Then, instead of having a probability of hopping along an edge, it has an amplitude of doing so. But as we’ve seen, a lot of similar math will still apply.

It should be fun to compare the two approaches. But I bet you’re wondering which approach is correct. This is a somewhat tricky question, at least for me. The answer would seem to depend on how much the ethyl cation is interacting with its environment—for example, bouncing off other molecules. When a system is interacting a lot with its environment, a probabilistic approach seems to be more appropriate. The relevant buzzword is ‘environmentally induced decoherence’.

However, there’s something much more basic I have tell you about.

After the paper by Balaban, Fǎrcaşiu and Bǎnicǎ came out, people gradually realized that the ethyl cation doesn’t really look like the drawing I showed you! It’s what chemists call ‘nonclassical’ ion. What they mean is this: its actual structure is not what you get by taking the traditional ball-and-stick model of an ethane molecule and ripping off a hydrogen. The ethyl cation really looks like this:

For more details, and pictures that you can actually rotate, see:

• Stephen Bacharach, Ethyl cation, Computational Organic Chemistry.

So, if we stubbornly insist on applying the Desargues graph to realistic chemistry, we need to find some other molecule to apply it to.

Trigonal bipyramidal molecules

Luckily, there are lots of options! They’re called trigonal bipyramidal molecules. They look like this:

The 5 balls on the outside are called ‘ligands’: they could be atoms or bunches of atoms. In chemistry, ‘ligand‘ just means something that’s stuck onto a central thing. For example, in phosphorus pentachloride the ligands are chlorine atoms, all attached to a central phosphorus atom:

It’s a colorless solid, but as you might expect, it’s pretty nasty stuff: it’s not flammable, but it reacts with water or heat to produce toxic chemicals like hydrogen chloride.

Another example is iron pentacarbonyl, where 5 carbon-oxygen ligands are attached to a central iron atom:

You can make this stuff by letting powdered iron react with carbon monoxide. It’s a straw-colored liquid with a pungent smell!

Whenever you’ve got a molecule of this shape, the ligands come in two kinds. There are the 2 ‘axial’ ones, and the 3 ‘equatorial’ ones:

And the molecule has 20 states… but only if count the states a certain way. We have to treat all 5 ligands as distinguishable, but think of two arrangements of them as the same if we can rotate one to get the other. The trigonal bipyramid has a rotational symmetry group with 6 elements. So, there are 5! / 6 = 20 states.

The transitions between states are devilishly tricky. They’re called pseudorotations, and they look like this:

If you look very carefully, you’ll see what’s going on. First the 2 axial ligands move towards each other to become equatorial.
Now the equatorial ones are no longer in the horizontal plane: they’re in the plane facing us! Then 2 of the 3 equatorial ones swing out to become axial. This fancy dance is called the Berry pseudorotation mechanism.

To get from one state to another this way, we have to pick 2 of the 3 equatorial ligands to swing out and become axial. There are 3 choices here. So, if we draw a graph with states as vertices and transitions as edges, it will have 20 vertices and 20 × 3 / 2 = 30 edges. That sounds suspiciously like the Desargues graph!

Puzzle 1. Show that the graph with states of a trigonal bipyramidal molecule as vertices and pseudorotations as edges is indeed the Desargues graph.

I think this fact was first noticed here:

• Paul C. Lauterbur and Fausto Ramirez, Pseudorotation in trigonal-bipyramidal molecules, J. Am. Chem. Soc. 90 (1968), 6722–6726.

Okay, enough for now! Next time I’ll say more about the Markov process or quantum process corresponding to a random walk on the Desargues graph. But since the Berry pseudorotation mechanism is so hard to visualize, I’ll pretend that the ethyl cation looks like this:

and I’ll use this picture to help us think about the Desargues graph.

That’s okay: everything we’ll figure out can easily be translated to apply to the real-world situation of a trigonal bipyramidal molecule. The virtue of math is that when two situations are ‘mathematically the same’, or ‘isomorphic’, we can talk about either one, and the results automatically apply to the other. This is true even if the one we talk about doesn’t actually exist in the real world!


Fool’s Gold

12 September, 2011

My favorite Platonic solids are the regular dodecahedron, with 12 faces and 20 corners:

and the regular icosahedron, with 12 corners and 20 faces:

They are close relatives, with all the same symmetries… but what excites me most is that they have 5-fold symmetry. It’s a theorem that no crystal can have 5-fold symmetry. So, we might wonder whether these shapes occur in nature… and if they don’t, how people dreamt up these shapes in the first place.

It’s widely believed that the Pythagoreans dreamt up the regular dodecahedron after seeing crystals of iron pyrite—the mineral also known as ‘fool’s gold’. Nobody has any proof of this. However, there were a lot of Pythagoreans in Sicily back around 500 BC, and also a lot of pyrite. And, it’s fairly common for pyrite to form crystals like this:


This crystal is a ‘pyritohedron’. It looks similar to regular dodecahedron—but it’s not! At the molecular level, iron pyrite has little crystal cells with cubic symmetry. But these cubes can form a pyritohedron:


(By the way, you can click on any of these pictures for more information.)

You’ll notice that the front face of this pyritohedron is like a staircase with steps that go up 2 cubes for each step forwards. In other words, it’s a staircase with slope 2. That’s the key to the math here! By definition, the pyritohedron has faces formed by planes at right angles to these 12 vectors:

\begin{array}{cccc} (2,1,0) &  (2,-1,0) &  (-2,1,0) &  (-2,-1,0) \\ (1,0,2) &  (-1,0,2) &  (1,0,-2) &  (-1,0,-2) \\ (0,2,1) &  (0,2,-1) &  (0,-2,1) &  (0,-2,-1) \\ \end{array}

On the other hand, a regular dodecahedron has faces formed by the planes at right angles to some very similar vectors, where the number 2 has been replaced by this number, called the golden ratio:

\displaystyle {\Phi = \frac{\sqrt{5} + 1}{2}}

Namely, these vectors:

\begin{array}{cccc} (\Phi,1,0) &  (\Phi,-1,0) &  (-\Phi,1,0) &  (-\Phi,-1,0) \\ (1,0,\Phi) &  (-1,0,\Phi) &  (1,0,-\Phi) &  (-1,0,-\Phi) \\ (0,\Phi,1) &  (0,\Phi,-1) &  (0,-\Phi,1) &  (0,-\Phi,-1) \\ \end{array}

Since

\Phi \approx 1.618...

the golden ratio is not terribly far from 2. So, the pyritohedron is a passable attempt at a regular dodecahedron. Perhaps it was even good enough to trick the Pythagoreans into inventing the real thing.

If so, we can say: fool’s gold made a fool’s golden ratio good enough to fool the Greeks.

At this point I can’t resist a digression. You get the Fibonacci numbers by starting with two 1’s and then generating a list of numbers where each is the sum of the previous two:

1, 1, 2, 3, 5, 8, 13, …

The ratios of consecutive Fibonacci numbers get closer and closer to the golden ratio. For example:

\begin{array}{ccl}  1/1 &=& 1 \\ 2/1 &=& 2  \\ 3/2 &=& 1.5   \\  5/3 &=& 1.6666..   \\ 8/5 &=& 1.6 \\ \end{array}

and so on. So, in theory, we could use these ratios to make cubical crystals that come closer and closer to a regular dodecahedron!

And in fact, pyrite doesn’t just form the 2/1 pyritohedron I showed you earlier. Sometimes it forms a 3/2 pyritohedron! This is noticeably better. The 2/1 version looks like this:

while the 3/2 version looks like this:

Has anyone ever seen a 5/3 pyritohedron? That would be even better. It would be quite hard to distinguish by eye from a true regular dodecahedron. Unfortunately, I don’t think iron pyrite forms such subtle crystals.

Okay. End of digression. But there’s another trick we can play!

These 12 vectors I mentioned:

\begin{array}{cccc} (\Phi,1,0) &  (\Phi,-1,0) &  (-\Phi,1,0) &  (-\Phi,-1,0) \\ (1,0,\Phi) &  (-1,0,\Phi) &  (1,0,-\Phi) &  (-1,0,-\Phi) \\ (0,\Phi,1) &  (0,\Phi,-1) &  (0,-\Phi,1) &  (0,-\Phi,-1) \\ \end{array}

besides being at right angles to the faces of the dodecahedron, are also the corners of the icosahedron!

And if we use the number 2 here instead of the number \Phi, we get the vertices of a so-called pseudoicosahedron. Again, this can be made out of cubes:



However, nobody seems to think the Greeks ever saw a crystal shaped like a pseudoicosahedron! The icosahedron is first mentioned in Book XIII of Euclid’s Elements, which speaks of:

the five so-called Platonic figures which, however, do not belong to Plato, three of the five being due to the Pythagoreans, namely the cube, the pyramid, and the dodecahedron, while the octahedron and the icosahedron are due to Theaetetus.

So, maybe Theaetetus discovered the icosahedron. Indeed, Benno Artmann has argued that this shape was the first mathematical object that was a pure creation of human thought, not inspired by anything people saw!

That idea is controversial. It leads to some fascinating puzzles, like: did the Scots make stone balls shaped like Platonic solids back in 2000 BC? For more on these puzzles, try this:

• John Baez, Who discovered the icosahedron?

But right now I want to head in another direction. It turns out iron pyrite can form a crystal shaped like a pseudoicosahedron! And as Johan Kjellman pointed out to me, one of these crystals was recently auctioned off… for only 47 dollars!

It’s beautiful:

So: did the Greeks ever seen one of these? Alas, we may never know.

For more on these ideas, see:

• John Baez, My favorite numbers: 5.

• John Baez, Tales of the dodecahedron: from Pythagoras to Plato to Poincaré.

• John Baez, This Week’s Finds in Mathematical Physics, "week241" and "week283".

• Ian O. Angell and Moreton Moore, Projections of cubic crystals, International Union of Crystallography.

To wrap up, I should admit that icosahedra and dodecahedra show up in many other places in nature—but probably too small for the ancient Greeks to see. Here are some sea creatures magnified 50 times:

And here’s a virus:


The gray bar on top is 10 nanometers long, while the bar on bottom is just 5 nanometers long.

The mathematics of viruses with 5-fold symmetry is fascinating. Just today, I learned of Reidun Twarock‘s recent discoveries in this area:

• Reidun Twarock, Mathematical virology: a novel approach to the structure and assembly of viruses, Phil. Trans. R. Soc. A 364 (2006), 3357-3373.

Most viruses with 5-fold symmetry have protein shells in patterns based on the same math as geodesic domes:

 

But some more unusual viruses, like polyomavirus and simian virus 40, use more sophisticated patterns made of two kinds of tiles:

They still have 5-fold symmetry, but these patterns are spherical versions of Penrose tilings! A Penrose tiling is a nonrepeating pattern, typically with approximate 5-fold symmetry, made out of two kinds of tiles:

To understand these more unusual viruses, Twarock needed to use some very clever math:

• Thomas Keef and Reidun Twarock, Affine extensions of the icosahedral group with applications to the three-dimensional organisation of simple viruses, J. Math. Biol. 59 (2009), 287-313.

But that’s another story for another day!



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 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.


Follow

Get every new post delivered to your Inbox.

Join 3,095 other followers