Chapter 1
PredatorPrey Dynamics for Rabbits, Trees, and Romance
Julien Clinton Sprott
Department of Physics
University of Wisconsin  Madison
sprott@juno.physics.wisc.edu
1.1. Introduction
The LotkaVolterra equations represent a simple nonlinear model for the dynamic interaction between two biological species in which one species (the predator) benefits at the expense of the other (the prey). With a change in signs, the same model can apply to two species that compete for resources or that symbiotically interact. However, the model is not structurally stable, since persistent timedependent (oscillatory) solutions occur for only a single value of the parameters.
This paper considers structurally stable variants of the LotkaVolterra equations with arbitrarily many species solved on a homogeneous twodimensional grid with coupling between neighboring cells. Interesting, biologicallyrealistic, spatiotemporal patterns are produced. These patterns emerge from random initial conditions and thus exhibit selforganization. The extent to which the patterns are selforganized critical (spatial and temporal scaleinvariant) and chaotic (positive Lyapunov exponent) will be examined.
The same equations, without the spatial interactions, can be used to model romantic relationships between individuals. Different romantic styles lead to different dynamics and ultimate fates. Love affairs involving more than two individuals can lead to chaos. Strange attractors resulting from such examples will be shown.
1.2. LotkaVolterra Equations
One variant of the LotkaVolterra equations (Murray 1993) for two species (such as rabbits and foxes) is
(1)
where R is the number of rabbits and F is the number of foxes, both positive and each normalized to its respective carrying capacity (the maximum allowed in the absence of the other), r_{1} and r_{2} are the respective growth rates in the absence of competition, and a_{1} and a_{2} determine the interspecies competition. In a predatorprey model, the predator (foxes) would have r_{2} < 0 and the other constants would be positive. However, the same equations with all positive constants could model competition, or with both growth rates negative could model cooperation or symbiosis (chickens and eggs, plants and seeds, bees and flowers, etc.).
1.3. Equilibrium and Stability
The system in Eq. (1) has four equilibria, one with no rabbits, one with no foxes, one with neither, and a coexisting one with
(2)
which is of primary interest.
For the predatorprey case (a_{1}r_{1} > 0, a_{1}r_{2} < 0), the coexisting equilibrium is a stable focus for r_{1}(1  a_{1}) < r_{2}(1  a_{2}), at which point it undergoes a Hopf bifurcation, after which the trajectory spirals outward without bound from the unstable focus. Hence there are no structurally stable oscillatory solutions. For the competition case (a_{1}r_{1} > 0, a_{1}r_{2} > 0), the coexisting equilibrium is a stable node for a_{1} < 1 and a_{2} < 1, at which point it undergoes a saddlenode bifurcation, after which one of the species dies while the other goes to its carrying capacity (R = 1 or F = 1). Thus there are no oscillatory solutions, and stability requires that the intraspecies competition dominates. When the interspecies competition dominates, the weaker species is extinguished by the principle of `competitive exclusion’ or `survival of the fittest’ (Gause 1971). For the cooperation case (a_{1}r_{1} < 0, a_{1}r_{2} < 0), both species either die or grow without bound.
With N species there are 2^{N} equilibria, only one of which represents coexistence, and it is unlikely that this equilibrium is stable since all of its eigenvalues must have have only negative real parts. Ecological systems exhibit diversity presumably because there are so many species from which to choose, because they are able to spread out over the landscape to minimize competition, and because species evolve to fill stable niches (Chesson 2000).. If many arbitrary species are introduced into a highly interacting environment, most would probably die.
1.4. Spatiotemporal Generalization
Now assume there are N species with population S_{i} for i = 1 to N and that they are spread out over a twodimensional landscape S_{i}(x,y). The species could be plants or animals or both. For convenience, take the landscape to be a square of size L with periodic boundary conditions, so that S_{i}(L,y) = S_{i}(0,y) and S_{i}(x,L) = S_{i}(x,0). One commonly assumes that each species obeys a reactiondiffusion equation of the form
(3)
This system is usually solved on a finite spatial grid of cells each of size d so that Ñ ^{2}S_{i}(x,y) = [S_{i}(x+d,y) + S_{i}(xd,y) + S_{i}(x,y+d) + S_{i}(x,yd) – 4S_{i}(x,y)] / d^{2}.
Diffusion is perhaps not the best model for biology, however, and if too large, it tends to produce spatially homogneity. Instead, assume that each species interacts not just with the other species in its own cell but also in the four nearestneighbor cells (a von Neumann neighborhood), giving
(4)
where = S_{j}(x+d,y) + S_{j}(xd,y) + S_{j}(x,y+d) + S_{j}(x,yd) + a _{j}S_{j}(x,y) is a weighted average of the neighborhood. In the example of rabbits and foxes, you can think of a _{j} as the tendency for the foxes to eat at home. With a _{j} = 0, the foxes always eat out, and with a _{j} = 1 they forage uniformly over a fivecell neighborhood. In the case of trees and seeds, this term is where one would include a seed dispersion kernel. The example that follows uses a _{j} = 1 for all j, but the results are not sensitive to the choice. Including only nearest neighbor cells normalizes space so that the cell size d is the order of the mean dispersal (or foraging) distance. Note that time can also be normalized to one of the growth times, so that we can take r_{1} = 1 without loss of generality.
1.5. Numerical Example
In Eq. (4) all the biology is contained in the vector r_{i}, the interaction matrix a_{ij}, and the dispersal vector a _{j} here taken as unitary. Instead of modeling realistic biology, we choose the values of r_{i} and a_{ij} from an IID random normal distribution with zero mean and unit variance and examine many instances of the model to explore a range of possible ecologies.
For brevity, Fig. 1 illustrates most of the common behaviors. It starts with six species with uniform random values S_{i}(x,y) in the range of 0 to 0.2 on a 100 ´ 100 grid. The upper plots show the spatial structure of the six species after 100 growth times, and the lower plot shows the cumulative relative abundance of each species versus time. One species (the fourth) dies out. The first, second, and sixth, nearly die, and then recover, after which the five species coexist with aperiodic temporal fluctuations and spatial heterogeneity.
Figure 1. Typical example of a spatiotemporal solution Eq. (4) with six initial species, one of which died.
Figure 2 shows the dominant species in each cell after 100 growth times, each in a different color. This display facilitates comparison with real data and with the results of cellular automata models (Sprott 2002). As in earlier studies, we define the cluster probability as the fraction of cells that are the same as their four nearest neighbors and Fourier analyze the temporal fluctuations in cluster probability to obtain its power spectrum. The result in Fig. 3 shows a power law over about a decade and a half, implying temporal scale invariance as suggestive of selforganized criticality (Bak 1996). Other quantities such as the total biomass S S_{i}(t) also have powerlaw spectra.
Figure 2. Landsccape pattern showing the dominant species in each cell.
Figure 3. Power spectrum of fluctuations in cluster probability.
To assess whether the dynamics are chaotic, we follow Lorenz (1963) and round the values of S_{i}(x,y) to four significant digits after the initial transient has decayed and calculate the growth of the error in the total biomass as the perturbed and unpertubed systems evolve deterministically. The result in Fig. 4 suggests an exponential growth in the error with a growth rate the order of 0.1r_{1}. If the system modeled a forest with a typical r_{1} of 50 years, the predictability time would be about 500 years. Five species appears to be the minimum number for such chaotic solutions. With four species, limit cycles were found, and with three or fewer species, all stable solutions appear to attract to a timeindependent equilibrium with no spatial structure. Spatial heterogeneity always correlates with temporal fluctuations.
Figure 4. Exponential growth in total biomass error suggesting chaos.
Note that the chaos and spatial structure arise from a purely deterministic model in which the only randomness is in the initial condition. In fact, similar structures arise from highly ordered initial conditions with noise as small as 10^{6}. The model is purely endogenous (no external effects), purely homogeneous (every cell is equivalent), and purely egalitarian (all species obey the same equation, with only different coefficients). The spatial patterns and fluctuations are inherent in the equations whose solutions spontaneously break the imposed symmetry.
To stress the generality of the model, we can apply it to romantic relationships. Imagine two lovers, Romeo and Juliet, characterized by a pair of equations such as Eq. (1) in which R is Romeo’s love for Juliet and F is Juliet’s love for Romeo, both positive. Each lover can be characterized by one of four romantic styles depending on the signs of r and a as shown in Table I using names adapted from Strogatz (1988). The variable r determines whether one’s love grows or dies in the absence of a response from the other (a = 0), and the variable a determines whether reciprocated love enhances or suppresses one’s feelings.
Table I. Romantic styles
a – + Cautious lover 
+ + Narcissicist nerd 
– – Hermit 
r + – Eager beaver 
With two interacting lovers, there are thus 2^{4} = 16 different combinations of romantic styles, the fate of which are determined by the strength of the interations. As an example, choosing r and a from an IID random normal distribution with zero mean and unit variance gives the results in Table II, where the percentages are the probability that a stable steady state is reached. In some sense, the best pairing is between an eager beaver and a narcissistic nerd, although two eager beavers have solutions that grow mutually without bound. Not surprisingly, the prospects are dismal for a hermit and not much better for a cautious lover. Fortunately, humans seem capable of adapting their romantic styles to fit the situation.
Table II. Probability of mutual stable love for various pairings.
Narcissistic nerd 
Eager beaver 
Cautious lover 
Hermit 

Narcissistic nerd 
46% 
67% 
5% 
0% 
Eager beaver 
67% 
39% 
0% 
0% 
Cautious lover 
5% 
0% 
0% 
0% 
Hermit 
0% 
0% 
0% 
0% 
It is also instructive to examine love triangles, in which case there are are four variables if two of the lovers are unaware of one another, and six variables if each person has feelings for the other two. The variables need not be romatic love ; the third person could be the child of a couple or perhaps a motherinlaw. With six variables, there are 2^{6} = 64 equilibria, only one of which represents a universally happy arrangement and 4^{6} = 4096 different combinations of styles, assuming each person can adopt a different interaction style toward each of the others. Not surprisingly, the prognosis for coexisting positive feelings is very low unless the individuals exhibit strong adaptability of their styles. Perhaps humans adapt to such situtions much the way plants and animals do, by limiting their interactions, evolving to fill a stable niche, drawing sustinence from others, and maintaining spatial separation.
With three or more variables, there is the possibility of chaotic solutions. However, a search for such solutions in the system of Eq. (1) generalized to six variables failed to reveal any such solutions, although other similar models do exhibit chaos (see for example http://sprott.physics.wisc.edu/lectures/love&hap/).
Acknowledgments
I would like to thank Warren Porter and Heike Lischke whose research inspired this work and George Rowlands and Janine Bolliger whose ongoing collaboration facilitated its completion.
References
Bak, P., 1996, How Nature Works: The Science of Selforganized Criticality, Corpernicus (New York).
Chesson, P., 2000, Mechanisms of Maintenance of Species Diversity. Annu. Rev. Ecol. Syst., 31, 343.
Gause, G. F., 1971, The Struggle for Existence, Dover (New York).
Lorenz, E.N., 1963, Deterministic Nonperiodic Flow. J. Atmos. Sci., 20, 130.
Murray, J. D., 1993, Mathematical Biology, Springer (Berlin).
Sprott, J. C., Bolliger, J., and Mladenoff, D. J., 2002, Selforganized Criticality in Forestlandscape Evolution. Phys. Lett. A, 297, 267.
Strogatz, S. H., 1988, Love Affairs and Differential Equations. Mathematics Magazine, 61, 35.