Chapter 1

The Lotka-Volterra 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 time-dependent (oscillatory) solutions occur for only a single value of the parameters.

This paper considers structurally stable variants of the Lotka-Volterra equations with arbitrarily many species solved on a homogeneous two-dimensional grid with coupling between neighboring cells. Interesting, biologically-realistic, spatio-temporal patterns are produced. These patterns emerge from random initial conditions and thus exhibit self-organization. The extent to which the patterns are self-organized critical (spatial and temporal scale-invariant) 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.

One variant of the Lotka-Volterra 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₁ and r₂ are the respective growth rates in the absence of competition, and a₁ and a₂ determine the interspecies competition. In a predator-prey model, the predator (foxes) would have r₂ < 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.).

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 predator-prey case (a₁r₁ > 0, a₁r₂ < 0), the coexisting equilibrium is a stable focus for r₁(1 - a₁) < -r₂(1 - a₂), 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₁r₁ > 0, a₁r₂ > 0), the coexisting equilibrium is a stable node for a₁ < 1 and a₂ < 1, at which point it undergoes a saddle-node 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₁r₁ < 0, a₁r₂ < 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.

Now assume there are N species with population S_i for i = 1 to N and that they are spread out over a two-dimensional 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 reaction-diffusion equation of the form

(3)

This system is usually solved on a finite spatial grid of cells each of size d so that Ñ ²S_i(x,y) = [S_i(x+d,y) + S_i(x-d,y) + S_i(x,y+d) + S_i(x,y-d) – 4S_i(x,y)] / d².

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 nearest-neighbor cells (a von Neumann neighborhood), giving

(4)

where = S_j(x+d,y) + S_j(x-d,y) + S_j(x,y+d) + S_j(x,y-d) + a _jS_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 five-cell 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 without loss of generality.

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 spatio-temporal 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 self-organized criticality (Bak 1996). Other quantities such as the total biomass S S_i(t) also have power-law 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₁. If the system modeled a forest with a typical r₁ 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 time-independent 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

With two interacting lovers, there are thus 2⁴ = 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.

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 mother-in-law. With six variables, there are 2⁶ = 64 equilibria, only one of which represents a universally happy arrangement and 4⁶ = 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 https://sprott.physics.wisc.edu/lectures/love&hap/).

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.