The Azimuth Project
Invasion fitness in moment-closure treatments (Rev #7, changes)

Showing changes from revision #6 to #7: Added | Removed | Changed

This page is a blog article in progress, written by Blake Stacey. To discuss as it is being written, go to the Azimuth Forum.


Moment closures are a way of forgetting information about a system in a controlled fashion, in the hope that an incomplete, fairly heavily “coarse-grained” picture of the system will still be useful in figuring out what will happen to it. Sometimes, this is a justifiable hope, but in other cases, we are right to wonder whether all the algebra it generates actually leads us to any insights. Here, we’ll be concerned with a particular application of this technology: studying the vulnerability of an ecosystem to invasion. We shall find expressions for invasion fitness, the expected relative growth rate of an initially-rare species or variety.

Chatty Presentation

Consider a lattice, each site of which can be occupied empty by an individual of “resident” type (RR), occupied by a mutant (MM), or empty (00 ), ). occupied The difference between the two is encoded in the choice of transition rules representing death, birth and migration. We can get an aggregate measure of the situation by finding the probability that a selfish randomly individual chosen ( site will be in stateSS) or occupied by an altruist (AA). The difference between selfishness and altruism is encoded in the choice of transition rules representing death, birth and migration. We can get an aggregate measure of the situation by finding the probability that a randomly chosen site will be in state aa, where aa can take values in the set {0R, A M,S0} \{0, \{R, A, M, S\} 0\}. A finer degree of distinction is provided by the conditional probabilities q a|bq_{a|b}, where, for example, q S R| A M q_{S|A} q_{R|M} denotes the probability that a randomly chosen neighbor site to a randomly chosen altruist mutant is selfish. of resident type. Note that if a selfish mutant is injected into a native resident population of altruists and its selfish offspring form a geographical cluster,q S M| S M q_{S|S} q_{M|M} can be much larger than p S M p_S p_M : few individuals are selfish mutants overall, but the probability of a selfish mutant life-form interacting with one another of mutant the same behavior type is high.

The pair dynamics of the system involves the time evolution of the probabilities p abp_{a b}, that is, the probability that a randomly selected lattice edge will have aa on one end and bb on the other. The differential equation for dp S R0M/dt d p_{S p_{R 0}/d M}/d t, for example, will have terms reflecting the processes which can form and destroy S R0M S0 R M pairs: 00RM S R0R 00\rightarrow R S M\rightarrow 0 R R , is one possibility, and S R0M00MM S R 0\rightarrow M\rightarrow 00 M M , is another. Death, which comes for organisms and leaves empty spaces behind, introduces processes like S R0M S RS0 S0\rightarrow R S M S \rightarrow R 0, S R S MS0M S R S\rightarrow M S \rightarrow 0 M , and S R0MS0A0 S R 0\rightarrow M S \rightarrow A 0 0 . and Reproduction can lead to formerly empty spaces becoming occupied: S RAS0RR S R A\rightarrow S 0 \rightarrow R R , let and us say. In language more like we’ve used in theM0MMM 0 \rightarrow M M. In language more like we’ve used in the Azimuth discussions of stochastic Petri net theory, we’ve moved beyond just creating and annihilating altruists residents and selfishists, mutants, and now we’re dynamically changing the number of “altruist–altruist” “resident–resident” and “altruist–selfish” “resident–mutant”pairs. So, we probably could make a Petri net diagram for these processes, but the labels on the boxes might start looking a little funny.

Each term in our differential equations will have a transition rate dependent upon a conditional probability of the form q a|bcq_{a|b c}, denoting the probability that a bb of a bcb c pair will have a neighbor of type aa. The differential equations for the pair probabilities p abp_{a b} thus depend on triplet probabilities p abcp_{a b c}, which depend upon quadruplet probabilities and so forth. To make progress, we truncate this hierarchy in a manner reminiscent to cutting off the BBGKY hierarchy in statistical mechanics: we impose the pair approximation that

q a|bcq a|b. q_{a|b c} \approx q_{a|b}.

This approximation destroys information about spatial structure and thereby introduces bias which in an ideal world ought to be accounted for.

(It just occurred to me that if a Petri net specifies a symmetric monoidal category, then successive truncations of the moment-dynamics hierarchy yield mappings between categories. Going from a pair approximation to a mean-field approximation, for example, transforms a Petri net whose circles are labelled with pair states to one labelled by site states. Category theory might be able to say something interesting here. Anything which can tame the horrible spew of equations which arises in these problems would be great to have.)

Invasion fitness is judged in the following manner. We start with a lattice devoid of altruists mutants (p A Ma=0 p_{A p_{M a} = 0) and find the equilibrium densities p S R S R * p_{S p_{R S}^* R}^* and p S R0 * p_{S p_{R 0}^* by setting

dp S R0dt=dp S R S Rdt=0. \frac{d p_{S0}}{d p_{R0}}{d t} = \frac{d p_{S p_{R S}}{d R}}{d t} = 0.

The exact form of p S R S R * p_{S p_{R S}^* R}^* and p S R0 * p_{S0}^* p_{R0}^* will depend upon interaction details which do we not won’t concern worry us about here. just yet. We then inject an a altruistic mutant strain into this situation; as the altruists mutants are initially rare, we can say they do not affect the large-scale dynamics of the selfish resident population. Summarizing the pair probabilitiesp A Ma p_{A p_{M a} with the shorthand p̲\underline{p}, we write the differential equation in matrix form

dp̲dt= M T(q a|bc)p̲, \frac{d\underline{p}}{d t} = M(q_{a|b T(q_{a|b c})\underline{p},

where the matrix M T(q a|bc) M(q_{a|b T(q_{a|b c}) encapsulates the details of our chosen behavior dynamics. modes. The pair approximation, in which we discard correlations of third and higher order, lets us simplify this to

dp̲dt= M T(q a|b)p̲. \frac{d\underline{p}}{d t} = M(q_{a|b})\underline{p}. T(q_{a|b})\underline{p}.

A When result people supported started by doing numerical simulations investigation of is lattice models like these, they found that the conditional probabilitiesq a| A M q_{a|A} q_{a|M} equilibrate. That is to say, even if the global density of altruists mutantsp A M p_A p_M changes, the local statistical structure of an a altruistic mutant cluster remains constant. This is the key statement which allows us to linearize the dynamics and write the behavior ofp̲\underline{p} in terms of eigenvectors and eigenvalues:

p̲(t)=Ce̲ Aexp(λt). \underline{p}(t) = C\underline{e}_A \exp(\lambda t).

The dominant eigenvalue λ\lambda of M T M T is the “invasion exponent” which characterizes whether an invasion will fail (λ<0\lambda \lt 0) or succeed (λ>0\lambda \gt 0). The eigenvector e̲ A\underline{e}_A associated with λ\lambda describes the vehicle of selection for the altruist’s mutants’ particular genetic variation, by summarizing the structure of the their altruists’ geographical cluster.

All of this, of course, is only as good as our linearization! If something interesting happens further away from the fixed point, just looking at the eigenvalues we got from our matrix MM won’t tell us about it.

All of this, of course, is only as good as our linearization! If something interesting happens further away from the fixed point, just looking at the eigenvalues we got from our matrix TT won’t tell us about it. In addition, if the actual dynamics of the system tend to form patterns which can’t be represented very well by pairwise correlations, then pair approximation will run into a wall. Minus van Baalen puts the issue in the following way:

The extent to which pair-dynamics models are satisfactory depends on the goal of the modeler. As we have seen, these models do not capture all of the phenomena that can be observed in simulations of fully spatial probabilistic cellular automata. Basically, the approximation fails whenever spatial structures arise that are difficult to “describe” using pairs alone. More technically, the method fails whenever significant higher-order correlations arise – that is, whenever the frequency of parituclar triplets (or triangles, squares, or all sorts of star-like configurations) starts to diverge from what one would expect on the basis of pair densities. Thus, pair-dynamics models satisfactorily describe probabilistic cellular automata in which only “small-scale” patterns arise. Larger, “meso-scale” patterns such as spirals are difficult to capture using this method.

—in Dieckmann et al. (2000), chapter 19.

A Specific Example 1: Birth, Death, Movement

I’ll probably take something from one of the essays in the Dieckmann et al. book. —Blake Stacey.

We follow van Baalen (in Dieckmann et al. (2000), chapter 19).

(Petri net pictures will go here.)

We write zz for the coordination number of the lattice. Birth:

R0RR, R 0 \rightarrow R R,

with rate b/zb/z.


Ra0a, R a \rightarrow 0 a,

with rate d/zd/z.

Movement or migration:

R00R, R 0 \rightarrow 0 R,

with rate m/zm/z.

dp R0dt=p R0[b/z+d+(11/z)q 0|R0m+(11/z)q R|0R(b+m)]+p 00(11/z)q R|00(b+m)+p RR[d+(11/z)q 0|RRm]. \frac{d p_{R 0}}{d t} = -p_{R 0}[b/z + d + (1 - 1/z)q_{0|R0}m + (1 - 1/z)q_{R|0R}(b + m)] + p_{00} (1 - 1/z) q_{R|00} (b + m) + p_{R R} [d + (1 - 1/z) q_{0|R R} m].
dp 00dt=p 002(11/z)q R|00(b+m)+p R02[d+(11/z)q 0|R0m]. \frac{d p_{0 0}}{d t} = -p_{00} 2(1 - 1/z) q_{R|00} (b + m) + p_{R 0} 2 [d + (1 - 1/z) q_{0|R 0} m].
dp RRdt=p R02[b/z+(11/z)q R|0R(b+m)]p RR2[d+(11/z)q 0|RRm]. \frac{d p_{R R}}{d t} = p_{R 0} 2[b/z + (1 - 1/z)q_{R|0 R} (b + m)] - p_{R R} 2[d + (1 - 1/z) q_{0|R R} m].
p i= jp ij. p_i = \sum_j p_{i j}.
dp Rdt=(bq 0|Rd)p R. \frac{d p_R}{d t} = (b q_{0|R} - d) p_R.

If we ignore spatial structure altogether, we can say that

q 0|R=p 0, q_{0|R} = p_0,

which by normalization of probability means

q 0|R=1p R. q_{0|R} = 1 - p_R.


dp Rdt=(b(1p R)d)p R. \frac{d p_R}{d t} = (b(1 - p_R) - d) p_R.

This should look familiar: it’s a logistic equation for population growth, with growth rate bdb - d and equilibrium population 1d/b1 - d/b.


  • H. Matsuda, N. Ogita, A. Sasaki, and K. Sato (1992), “Statistical mechanics of population”, Progress of Theoretical Physics 88, 6: 1035–49. (web)
  • U. Dieckmann, R. Law, and J. A. J. Metz, eds., The Geometry of Ecological Interactions: Simplifying Spatial Complexity. Cambridge University Press, 2000.
  • B. Allen, Studies in the Mathematics of Evolution and Biodiversity. PhD thesis, Boston University, 2010. (web)
  • J. A. Damore and J. Gore (2011), “Understanding microbial cooperation”. Journal of Theoretical Biology DOI:10.1016/j.jtbi.2011.03.008 (pdf).
  • B. Simon, J. A. Fletcher and M. Doebeli (2011), “Hamilton’s rule in multi-level selection models” Journal of Theoretical Biology (web)
  • B. C. Stacey, A. Gros and Y. Bar-Yam (2011), “Beyond the mean field in host-pathogen spatial ecology”, arXiv:1110.3845 [nlin.AO].

category: blog