Equilibrium and stability are very important concepts in ecology, but they have many conflicting definitions. One of the most used was brought from the branch of physics and mathematics called dynamical systems analysis.

It is this approach that brought to ecology equations to describe the dynamics of populations, such as logistic growth equation and the system of Lotka-Volterra equations.

There are techniques designed to evaluate whether these systems of equations have equilibrium points, and, if so, whether these equilibriums are stable. Robert May (1972) used these tools to demonstrate a somewhat surprising result: the stability in food webs decreases with the increase of species and complexity.

This exercise is a simplified repetition of the procedure used by Robert May. The goal is that you understand the concepts of equilibrium and stability in dynamical systems, and to differentiate them from other equilibrium and stability concepts used in ecology. With this you will also have the elements to critically evaluate the results of May (1972) work.

This exercise will be done using R (R Core Team 2012), but you don't need to know the R language, because we will provide the commands ready to be executed. They are printed in this page, but also on the annex script. The only thing you need to know is how to run these commands on R. You can either copy and paste them in the R command line, or use the script file. Just follow these steps:

1. Install the R environment in your computer, along with the libraries deSolve and rootSolve. The offical R page has some instructions, but you can also read our

R installation tutorial.

1. Create a new folder to place your exercises.
2. Copy to this folder the following files:
   1. eq_funcoes.r
   2. eq_LV_comandos.r
3. Open R and the script file eq_comandos.r. Make sure that your working directory is where it the files are.
4. The commands in the script are in the same order as they are presented in this tutorial. Follow the tutorial, sending the commands as they are required.
5. If you're not sure how to send the commands, follow this tutorial.
6. Load in your R session the packages and functions that we will use in this tutorial:

library(deSolve)
library(rootSolve)
source("eq_funcoes.r")

In the stability tutorial, we have seen the criteria for local stability in dynamic systems:

The first step to understand Robert May's results is to generalize this criterion for more interacting populations. First, let's use Lotka-Volterra's competitive equation system:

$$V_1 \ = \ r_1 N_1 \left( \frac{K_1 - N_1 - \alpha N_2}{K_1} \right)$$ $$V_2 \ = \ r_2 N_2 \left( \frac{K_2 - N_2 - \beta N_1}{K_2} \right)$$

Here, the variables represent, for species 1 and 2:

• $V_1$, $V_2$ are the population growth velocities
• $r_1$, $r_2$ are the intrinsic growth rates
• $K_1$, $K_2$ are the support capacities
• $\alpha$, $\beta$ are the competition coefficients of each species over the other

Use the function plota.LV to plot the abundances of both competing species. The arguments for this function are the parameters explained above, plus the initial condition:

• n : vector with the initial population sizes for both species
• time : maximum simulation time

Verify that the function is working with the following command:

## First combination of parameters: 1/beta < k1/k2 < alfa  ##
plota.LV(n=c(n1=1,n2=1),r1=0.2,K1=150,r2=0.2,K2=100,
  alfa=0.2,beta=0.1,time=150)

In this system, the equilibrium population sizes are (Gotelli 2007):

$$\hat N_1 = \frac{K_1-\alpha K_2}{1-\alpha \beta}$$
$$\hat N_2 = \frac{K_2-\beta K_1}{1-\alpha \beta}$$

The function lv.neq calculates this values. Use it to obtain the equilibrium values for the same parameter set, and compare the value with the graph you have produced:

(lv.n1 <- lv.neq(r1=0.2,K1=150,r2=0.2,K2=100,alfa=0.2,beta=0.1))

There are two conditions for the coexistence of the competing species in Lotka-Volterra's system:

$$\frac{1}{\beta} \ > \ \frac{K_1}{K_2} \ > \ \alpha$$
$$\frac{1}{\beta} \ < \ \frac{K_1}{K_2} \ < \ \alpha$$

Let's see if the equilibrium population sizes are stable. The function plota.LV has also parameters to introduce perturbations in a specified time.

We have already calculated the equilibrium point for our parameter set, and stored its values in the object lv.n1. To make the graph now we can use this values as the parameter n of the function. But now, use the argument perturb to introduce a perturbation, subtracting one individual from species 1 and introducing one for species 2:

## Perturbation from equilibrium ##
plota.LV(n=lv.n1,r1=0.2,K1=150,r2=0.2,K2=100,alfa=0.2,beta=0.1,time=150, perturb=c(-1,1), t.perturb=50)

Do the same for a combination of parameters that satisfies the second condition above:

## Coexistence of two species, 1/beta > k1/k2 > alfa ##
## Equilibrium populations
(lv.n2 <- lv.neq(r1=0.2,K1=150,r2=0.2,K2=100,alfa=1.8,beta=0.9))
## Graph ##
plota.LV(n=lv.n2,r1=0.2,K1=150,r2=0.2,K2=100,alfa=1.8,beta=0.9,time=150)

## Perturbation from equilibrium ##
plota.LV(n=lv.n2,r1=0.2,K1=150,r2=0.2,K2=100,alfa=1.8,beta=0.9,time=200, perturb=c(-1,1), t.perturb=50)

What happens if the perturbation direction is reversed?

## Reversing the perturbation ##
plota.LV(n=lv.n2,r1=0.2,K1=150,r2=0.2,K2=100,alfa=1.8,beta=0.9,time=200, perturb=c(1,-1), t.perturb=50)

The logic behind the local stability of this system is similar to the stability of the model with one population. One point is locally stable if, on the neighboorhood of this point, the speed (variation of the population size) has a negative relation to the population size. But now we have two equations, and in each of them the velocities depend on both population sizes.

One way to solve this question is by evaluating the effect of each population over the velocities, keeping the other as fixed. This is done using partial derivatives, which are represented by the symbol $\partial$. For example:

$$\frac{\partial V_1}{\partial N_1}$$

This is the derivative of the growth velocity of population 1 in respect to the size of population 1, and keeping population 2 constant. This derivative represents the effect that changing this population size has over its own population growth. The effect that population 1 has over population 2 is given by:

$$\frac{\partial V_2}{\partial N_1}$$

Next trick we can use is organizing this partial derivatives in a matrix, that we can interpret as having the combinations of the effects that every species has on itself and on the competitor:

$$\left( \begin{matrix} \frac{\partial V_1}{\partial N_1} & \frac{\partial V_1}{\partial N_2} \\ & \\\frac{\partial V_2}{\partial N_1} & \frac{\partial V_2}{\partial N_2} \end{matrix} \right) $$

This is called the Jacobian matrix of the equation system. The diagonal of the matrix contains the effects that every population has on itself, and off-diagonals are the inter-specific terms.

Studying the local stability of the logistic equation, we have evaluated the sign of the derivative in the equilibrium points. Now, we do the same thing for each partial derivative. Doing so, we will end up having the matrix of the partial derivatives, evaluated on the equilibrium point. May designated this Jacobian matrix for the Lotka-Volterra system evaluated on the equilibrium point as the Community matrix.

For one equation, the linear approximation indicates that the equilibrium is stable if the relation between the velocities and population sizes in its neighborhood is negative. But how can we assess that now? The extension for the stability criterion for a Lotka-Volterra system is as follows:

Eigenvalues are a property of every matrix, and their study is one of the basis of linear algebra. For a didactic introduction, see Otto & Day (2007, caps 7 e 8).

For our purposes, it needs only be stated that negative eigenvalues are the multivariate correspondent of the negative derivative for one variable.

The number of eigenvalues for a community matrix is equal to the number of species. The calculation of eigenvalues is not a simple task, but we have R to help us with its eigen function.

Let's calculate the community matrix for our first parameter set. To do this, use the jacob.lv function, informing the value of the parameters. This function finds out the equilibrium point and calculates the Jacobian matrix on those points:

## First combination of parameters: 1/beta < k1/k2 < alfa  ##
## Community matrix
j1 <- jacob.lv(r1=0.2,K1=150,r2=0.2,K2=100,alfa=0.2,beta=0.1)

Now we apply the function eigen to obtain the eigenvalues of this matrix:

##Eigenvalues
eigen(j1, only.values=TRUE)$values

Do the same for the second parameter combination we had:

## Second parameter combination: 1/beta > k1/k2 > alfa  ## 
## Community matrix
j2 <- jacob.lv(r1=0.2,K1=150,r2=0.2,K2=100,alfa=1.8,beta=0.9)
##Eigenvalues
eigen(j2, only.values=TRUE)$values

Question: you have already introduced a perturbation in these systems. Are the results for the community matrix eigenvalues coherent with what you observed after the perturbations?

The community matrix can be extended to how many species we wish. Robert May used this to assess whether an increase in species diversity increases stability. If this is true, as was accepted at the time, it should be easier to assemble stable matrices that are large than small.

To test this hypothesis, May made the following simulation:

1. Determine the number of system species, S
2. Determine the fraction of potential interactions that occur. This is the connectance matrix, C
3. Create an array with S lines and S columns
4. Complete the diagonal of this matrix with -1
5. Fill at random the remaining cells with random values until it reaches the desired connectance
6. Check whether the eigenvalues of this matrix are all negative. If so, count this as a stable system

Thus, we generate a system in which there is the same amount of intra-specific competition for all species and a fixed proportion of inter-specific associations, distributed at random. The value and sign of these interactions is also drawn at random. May randomly selected inter-specific effects of a zero mean normal distribution, so that positive and negative interactions are equally likely. Another important parameter of the simulation is the standard deviation from this normal distribution, which increases the average value of the interaction effect. Therefore, it is called interaction force.

Our function may performs a simulation according to the rules above. Its parameters are:

• S : number of species
• C : connectance
• f : mean interaction force
• nsim: number of repetitions

The function returns the proportion of repetitions that produced stable systems¹⁾.

Let's start with 20 species, connectance 0.3 ²⁾ and 0.2 interaction force:

(sim.1 <- may(S=20,C=0.3, f=0.2, nsim=100))

The command above stores the results in a table. Continue by increasing the number of species, and keeping the other parameters constant:

(sim.1 <- rbind(sim.1,may(S=40,C=0.3, f=0.2, nsim=100)))
(sim.1 <- rbind(sim.1,may(S=60,C=0.3, f=0.2, nsim=100)))
(sim.1 <- rbind(sim.1,may(S=80,C=0.3, f=0.2, nsim=100)))
(sim.1 <- rbind(sim.1,may(S=90,C=0.3, f=0.2, nsim=100)))
(sim.1 <- rbind(sim.1,may(S=100,C=0.3, f=0.2, nsim=100)))
(sim.1 <- rbind(sim.1,may(S=110,C=0.3, f=0.2, nsim=100)))
(sim.1 <- rbind(sim.1,may(S=120,C=0.3, f=0.2, nsim=100)))

What can you see here? A graph may help:

## Graph
plot(p.estab~S, data=sim.1, xlab="N of species",ylab="Stable matrices proportion")

The connectance can be interpreted as a measure of system complexity. Does greater complexity increases the chances of stability? Try it out:

## Increasing the connectance ##
(sim.2 <- may(S=120,C=0.3, f=0.2))
(sim.2 <- rbind(sim.2,may(S=120,C=0.28, f=0.2)))
(sim.2 <- rbind(sim.2,may(S=120,C=0.26, f=0.2)))
(sim.2 <- rbind(sim.2,may(S=120,C=0.24, f=0.2)))
(sim.2 <- rbind(sim.2,may(S=120,C=0.22, f=0.2)))
(sim.2 <- rbind(sim.2,may(S=120,C=0.20, f=0.2)))
(sim.2 <- rbind(sim.2,may(S=120,C=0.18, f=0.2)))
(sim.2 <- rbind(sim.2,may(S=120,C=0.16, f=0.2)))
## Graph
plot(p.estab~C, data=sim.2, xlab="Connectance",ylab="Stable matrices proportion")

Question: do diversity or complexity lead to an increase in stability?

• Gotelli, N. 2007. Ecologia. Londrina, Ed. Planta. (A basic reference about the dynamic models in ecology).
• May, R.M. 1972. Will a large complex system be stable? Nature, 238, 413-414. (The classic paper that established the concept of trophic network stability as a solution for a Lotka-Volterra system.)
• May, R.M. 2001. Stability and complexity in model ecosystems. Princeton, Princeton University Press. (In this influent monograph, Robert May develops the ideas of the 1972 paper. The first edition is from 1973, but it was republished by Princeton Landmarks of Biology in 2001.)
• Sarah P. Otto & Troy Day 2007. A Biologist's Guide to Mathematical Modeling in Ecology and Evolution. Princeton, Princeton University Press. (A great introduction to the maths behind these models, from biologists to biologists, using more intuitive approaches to the problems. A great source if you want to understand the details of the analyses used in this tutorial. See also the book site.)
• R Development Core Team (2012). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org/.