spatially explicit metapopulation models in scipy.

Making Pretty Pictures

I started to learn to program about 5 years ago running population ecology models in mathematica. Yesterday, I found an old mma notebook with a model modified to include differential parasitoid dispersal to adjacent host cells depending on the host density in those cells. It's bascially a discrete-time Nicholson-Bailey model. But in a grid of cells, where each cell contains a population of hosts (H) and parasitoids (P) that give birth, die, eat, and get eaten according to the NB model. Each generation, following birth/reproduction/predation, the hosts and parasitoids disperse. The hosts disperse equally to the 8 surrounding cells in their neigbhorhood. The parasitoids can move irrespective of host densities when the aggregation parameter (eta) is 0. When aggregation is 1 (eta == 1), the parasitoids move to adjacent cells in exact proportion to host densities in each of the surrounding cells. muH and muP (i'm too lazy to figure out how to write the symbols for mu) determine the proportion of individuals in each population that disperse. I no longer have mathematica, but I downloaded all 96 megabytes of the reader to open the model:

I'm no longer a big fan of mathematica, but, it does make it easy to use funky characters and show equations nicely.
There, the ListCorrelate does the dispersal. It's a convolution to "smear" out populations according to the kernels, which are defined by the muH, muP and for the parasitoid, eta.
It took literally 10 minutes to translate that to numpy/scipy.
I'll paste that code at the end.
The cool thing is that you can see complex spatial patterns arising within the gridded metapopulation, even when the sum of the densities across cells is constant. So, running the model, and plotting the time series (top), and the phase portrait (bottom) are shown.

The time-series is the average density of Hosts, Parasitoids across the metaopulation. It takes a while to stabilize, but they oscillate down to an equilibrium. The phase portrait is just plotting the mean density of hosts vs parasitoids for each generation. The generations are discrete, so the lines are merely to view the trajectory. It's hard to see from the scale of the time series, but the parasitoid cycles lag those of the host by about a generation.
But the model is running in a spatially explicit manner, and averaging loses all of that spatial data. So, starting at 200 generations, the program then saves a snapshot of the grid, with higher host densities in white, and lower values in black. This is done every third generation. After the model run, there's a directory with a bunch of images. It'd be nice to play them as a movie... Imagemagick's convert is perfect for this:

convert -delay 40 -loop 0 images/*.png metapop.gif

Where that command creates a nice little (4MB) animated gif:

That shows the spatiotemporal dynamics of the host population.
That's a pretty simple way to incorporate time into an visualization when you're otherwise limited to 2 dimensions... Maybe I'm too easily amused, but I could watch that all day.
Anyway, here's the hastily translated code:

import numpy
import pylab
from scipy.signal import convolve2d

def doplot(hm, pm):
    pylab.plot(range(gens), hm, 'b')
    pylab.plot(range(gens), pm, 'g')
    pylab.ylim(0, pm.mean() + hm.mean())
    pylab.legend(('hosts', 'parasitoids'))
    pylab.subplot(212)
    pylab.plot(hm, pm, 'r')
    pylab.plot(hm, pm, 'k.')
    pylab.xlim(0, 1.3 * numpy.mean(hm) + 5)
    pylab.ylim(0, 1.3 * numpy.mean(pm) + 5)
    pylab.xlabel('hosts')
    pylab.ylabel('parasitoids')
    pylab.show()

def metapop(H0, P0, a=0.05, l=3., K=1000, muH=0.2, muP=0.8
              , Hrange=1, Prange=1, eta=3, size=32, gens=1000
              , mode="same", boundary="wrap"):
    E = numpy.e
    Hmeta = []
    Pmeta = []
    Pkern = numpy.ones((2 * Prange + 1, 2 * Prange + 1)
            , dtype=numpy.double)
    Hkern = numpy.ones((2 * Hrange + 1, 2 * Hrange + 1)
            , dtype=numpy.double)
    Hkern *= muH/((2 * Hrange + 1)**2 - 1.)
    Hkern[Hrange + 1, Hrange + 1] = 1.0 - muH

    for gen in range(1, gens + 1):
        Hmeta.append(H0.mean())
        Pmeta.append(P0.mean())
        
        # poisson search.
        f = E**(-a * P0)

        # predation, births
        H1 = l * H0 * f * numpy.exp(-numpy.log(l) * H0*f/K)
        P1 = H0 * 1 - f

        # simple movement between adjacent cells.
        H0 = 0.001 + convolve2d(H1, Hkern, mode=mode, boundary=boundary)

        # biased movement by parasitoid according to density of hosts in
        # adjacent cells.
        heta = H0**eta
        B = heta / convolve2d(heta, Pkern, mode=mode
               ,  boundary=boundary)

        P0 = (1. - muP) * P1 + muP * B * convolve2d(P1, Pkern, mode=mode
                , boundary=boundary)
        P0 *= P1.sum()/P0.sum()
        if gen > 200 and not gen % 3:
            pylab.figure(figsize=(4,4))
            pylab.axes([0, 0, 1, 1])
            pylab.imshow(H0, cmap=pylab.cm.gray)
            pylab.xticks([])
            pylab.yticks([])
            pylab.savefig('images/%03i.png' % gen)
            pylab.close()

    return numpy.array(Hmeta), numpy.array(Pmeta)

if __name__ == "__main__":
    gens = 300
    size = 64
    H0 = 1000 * numpy.abs(numpy.random.randn(size * size).reshape(size, size))
    P0 = numpy.zeros_like(H0)
    P0[0, 0] = 1.
    P0[32, 32] = 1.

    pylab.close()
    pylab.subplot(211)
    hm, pm = metapop(H0, P0, gens=gens, a = 0.02)
    doplot(hm, pm)

The variable names are ugly, but it's a start for something to play with. It's interesting to see the effects of parasitoid aggregation (eta). And the parasitoid initialization--here, it introduces a single parasitoid at cells 0,0 and 32, 32. Should anyone actually use this, I have a derivation of the equations, and better explanation of the parameters available.