✎ Diversity and Stability []

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
<WRAP tabs>
  * [[en:ecovirt:roteiro:sucess:div_estab|{{:ecovirt:rlogo.png?20|}}]]
</WRAP>

====== Diversity and Stability ======
{{http://imgs.xkcd.com/comics/computer_problems.png?400}}

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 [[http://en.wikipedia.org/wiki/Dynamical_system|dynamical systems analysis]].

It is this approach that brought to ecology equations to describe the dynamics of populations, 
such as [[http://en.wikipedia.org/wiki/Logistic_function#In_ecology:_modeling_population_growth|logistic growth equation]] and the system of [[http://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equation|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.

<WRAP center round important 60%>
In order to follow this tutorial, it is assumed that you are familiar with the criteria for analysis of equilibria and their stability in population models.
If that is not the case, first follow the tutorial on [[en:ecovirt:roteiro:math:stabilitysage|stability in dynamical systems]]
</WRAP>

===== Preparation: the R environment =====
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:

  - Install the R environment in your computer, along with the libraries ''deSolve'' and ''rootSolve''. The [[http://www.r-project.org/|offical R page]] has some instructions, but you can also read our 
[[en:ecovirt:roteiro:soft:instalacaor| R installation tutorial]].
  - Create a new folder to place your exercises.
  - Copy to this folder the following files:
    - {{:ecovirt:roteiro:math:eq_funcoes.r}}
    - {{:ecovirt:roteiro:sucess:eq_lv_comandos.r|}}
  - Open R and the script file ''eq_comandos.r''. Make sure that your working directory is where it the files are.
  - 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.
  - If you're not sure how to send the commands, follow this [[en:ecovirt:roteiro:soft:rprincip|tutorial]].
  - Load in your R session the packages and functions that we will use in this tutorial:

<code>
library(deSolve)
library(rootSolve)
source("eq_funcoes.r")
</code>

===== Two populations =====

In the [[en:ecovirt:roteiro:math:stabilitysage|stability tutorial]], we have seen the criteria for local stability in dynamic systems:

<WRAP center round box 60%>
One equilibrium point for a population dynamics system is stable if the derivative of the population growth rate is negative in that point.
</WRAP>

The first step to understand Robert May's results is to generalize this criterion for more interacting populations.
First, let's use [[http://en.wikipedia.org/wiki/Competitive_Lotka%E2%80%93Volterra_equations|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:

<code>
## 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)
</code>

==== Equilibrium and Coexistence ====

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:

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

==== Stability ====

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:

<code>
## 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)
</code>

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

<code>
## 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)
</code>

What happens if the perturbation direction is reversed?

<code>
## 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)
</code>

==== Mathematical interpretation====

The logic behind the local stability of this system is similar to the stability of the model with [[en:ecovirt:roteiro:math:stabilitysage#estabilidade_da_logistica|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 [[http://en.wikipedia.org/wiki/Partial_derivative|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:

<WRAP box red 50% center round>
For a Lotka-Volterra system of equations, one equilibrium point is stable if all the eigenvalues of the community matrix are negative.
</WRAP>

[[http://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors|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 [[http://stat.ethz.ch/R-manual/R-patched/library/base/html/eigen.html|eigen]] function.

=== R Calculations ===

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:

<code>
## 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)
</code>

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

<code>
##Eigenvalues
eigen(j1, only.values=TRUE)$values
</code>

Do the same for the second parameter combination we had:

<code>
## 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
</code>

**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?

===== Baron May's Work =====

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:
  - Determine the number of system species, S
  - Determine the fraction of potential interactions that occur. This is the connectance matrix, C
  - Create an array with S lines and S columns
  - Complete the diagonal of this matrix with -1
  - Fill at random the remaining cells with random values until it reaches the desired connectance
  - 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**.

==== Simulating ====

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((i.e., arrays of communities in which the eigenvalues are all negative)).

Let's start with 20 species, connectance 0.3 ((30% of the community matrix cells have non-zero value)) and 0.2 interaction force:

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

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

<code>
(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)))
</code>

What can you see here? A graph may help:
<code>
## Graph
plot(p.estab~S, data=sim.1, xlab="N of species",ylab="Stable matrices proportion")
</code>

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

<code>
## 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")
</code>

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

===== To learn more =====

  * **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 [[http://press.princeton.edu/catalogs/series/plb.html|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 [[http://www.zoology.ubc.ca/biomath/|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/]].
  
{{tag> R cálculo derivada equação_diferencial lotka-volterra crescimento_logístico}}
User Tools

Site Tools

Page Tools
