Stability and complexity in model ecosystems Level 2 module in "Modelling course in population and evolutionary biology" (701-1418-00SS) Module author:Sebastian Bonhoeffer Course director:Sebastian Bonhoeffer Theoretical Biology Institute of Integrative Biology ETH Zuirich 1 Introduction The title of this module is in fact the title of an influential book that Robert M.Maya wrote in the 1970'sb.In this book May addressed the relationship between community structure in ecosystems and their properties of stability.At the time considerable evidence had been accumulated to suggest that in nature structurally complex communities tend to be more stable towards perturbations than simple ones.These observations led to an almost dogmatic view that stability results from the complexity of natural ecosystems.May challenged this view in the 1970'se by showing mathematically that in model ecosystems stability decreases with the number of species in a community and can therefore not be regarded as a straightforward consequence of complexity.To derive this result May investigated the stability properties of *Robert M.May (now Lord May of Oxford)has had an outstanding impact on theoretical ecology.He has made seminal contributions to many areas such as chaos in biology (see also the module on the logistic difference equation),theoretical epidemiology (see modules on infectious diseases),and to stability of ecosystems (this module).Being a native of Australia,Robert May was first a theoretical physicist.Becoming the successor to Robert McArthur as a theoretical ecologist at Princeton University,Robert May moved on later to Oxford University.As one of very few scientist Robert May was made a Lord by the British Queen in 2001.Robert May has also received an honorary doctorate from the ETH Zurich. PR.M.May,"Stability and complexity in model ecosystems",2nd ed,Princeton University Press,1974 In addition to the book mentioned above,May published an article in Nature in 1972 entitled "Will a large complex system be stable?
Stability and complexity in model ecosystems Level 2 module in “Modelling course in population and evolutionary biology” (701-1418-00 SS) Module author: Sebastian Bonhoeffer Course director: Sebastian Bonhoeffer Theoretical Biology Institute of Integrative Biology ETH Z¨urich 1 Introduction The title of this module is in fact the title of an influential book that Robert M. Maya wrote in the 1970’sb . In this book May addressed the relationship between community structure in ecosystems and their properties of stability. At the time considerable evidence had been accumulated to suggest that in nature structurally complex communities tend to be more stable towards perturbations than simple ones. These observations led to an almost dogmatic view that stability results from the complexity of natural ecosystems. May challenged this view in the 1970’sc by showing mathematically that in model ecosystems stability decreases with the number of species in a community and can therefore not be regarded as a straightforward consequence of complexity. To derive this result May investigated the stability properties of aRobert M. May (now Lord May of Oxford) has had an outstanding impact on theoretical ecology. He has made seminal contributions to many areas such as chaos in biology (see also the module on the logistic difference equation), theoretical epidemiology (see modules on infectious diseases), and to stability of ecosystems (this module). Being a native of Australia, Robert May was first a theoretical physicist. Becoming the successor to Robert McArthur as a theoretical ecologist at Princeton University, Robert May moved on later to Oxford University. As one of very few scientist Robert May was made a Lord by the British Queen in 2001. Robert May has also received an honorary doctorate from the ETH Zurich. bR.M. May, “Stability and complexity in model ecosystems”, 2nd ed, Princeton University Press, 1974 c In addition to the book mentioned above, May published an article in Nature in 1972 entitled “Will a large complex system be stable?”. 1
randomly assembled interacting species by means of eigenvalue analysis of random interaction matrices.Specifically he showed that the likelihood that an ecosystem is stable decreased with the fraction of interactions between species that are realized (connectivity)and the strength of these interactions. May's work opened the door towards a more formal investigation of the relationship between ecosystem structure and stability.However,May's approach also has some limitations.May studied randomly assembled communities,while real communities are likely not random,and neither are species extinctions:clearly,rare species are more likely to become extinct than abundant ones.Until today the relationship between community structure and stability remains one of the major unresolved questions in ecology. 2 Multi-species Lotka-Volterra models In this module we do not follow May's original approach,but instead will numerically simulate multi-species Lotka-Volterra models to investigate how ecosystem stability is related to size. The Lotka-Volterra model (LVM)for two speciesd is defined by dm/dt=r[1-(n1+12n2)/K1]n1 (1) dn2/dt=r2[1-(n2+21n1)/K2]n2 (2) where ni and n2 are the densities of populations 1 and 2.The interpretation of the parameters r1.2 and K1.2 are equivalent to the logistic differential equation,i.e.r1.2 are the maximal growth rates and K1.2 are the carrying capacities for the respective species.The interpretation of parameters Y12 and 721 becomes clear from inspecting the above equations.They describe the effects of competition between the two populations.If 2 >1,this implies that the negative effect of species 2 on species 1 is stronger than the negative effect of species 1 on itself.Hence, in this case for species 1 intraspecific competition is weaker than interspecific competition. Conversely,if y120 then the species j inhibits the growth of species i.Furthermore,if dSee also the script ofS.Bonhoeffer's lecture"Ecology and Evolution II:Populations"on the main course web page. 2
randomly assembled interacting species by means of eigenvalue analysis of random interaction matrices. Specifically he showed that the likelihood that an ecosystem is stable decreased with the fraction of interactions between species that are realized (connectivity) and the strength of these interactions. May’s work opened the door towards a more formal investigation of the relationship between ecosystem structure and stability. However, May’s approach also has some limitations. May studied randomly assembled communities, while real communities are likely not random, and neither are species extinctions: clearly, rare species are more likely to become extinct than abundant ones. Until today the relationship between community structure and stability remains one of the major unresolved questions in ecology. 2 Multi-species Lotka-Volterra models In this module we do not follow May’s original approach, but instead will numerically simulate multi-species Lotka-Volterra models to investigate how ecosystem stability is related to size. The Lotka-Volterra model (LVM) for two speciesd is defined by dn1/dt = r1[1 − (n1 + γ12n2)/K1]n1 (1) dn2/dt = r2[1 − (n2 + γ21n1)/K2]n2 (2) where n1 and n2 are the densities of populations 1 and 2. The interpretation of the parameters r1,2 and K1,2 are equivalent to the logistic differential equation, i.e. r1,2 are the maximal growth rates and K1,2 are the carrying capacities for the respective species. The interpretation of parameters γ12 and γ21 becomes clear from inspecting the above equations. They describe the effects of competition between the two populations. If γ12 > 1, this implies that the negative effect of species 2 on species 1 is stronger than the negative effect of species 1 on itself. Hence, in this case for species 1 intraspecific competition is weaker than interspecific competition. Conversely, if γ12 0 then the species j inhibits the growth of species i. Furthermore, if dSee also the script of S. Bonhoeffer’s lecture ”Ecology and Evolution II: Populations” on the main course web page. 2
aj>1 then the effect of species j on i is stronger that the effect of species i on itself.If all aij >0 then the multispecies LVM describes a community with only competitive interactions (i.e.without predator-prey relations between species).Predator-prey relations are characterized by interactions where aij0.(Here species i is a predator of species j).Note, however,that for reasons of simplicity we will focus only on competitive interactions.(One of the advanced topics suggested below is to extend the approach to include predator-prey type interactions.) 2.1 Numerical integration of the multi-species LVM in R The accompanying starting R script provides an implementation for the numerical solution of the multi-species LVM.Please consult this script in parallel to reading this section.The R script is commented extensively.The numerical simulation of the multi-species LVM can be programmed in a few lines in R(if you know how!).Essentially what you need is a function that defines the derivative of the LVM(i.e.the right hand side of equation 3).This function looks like this: ### #lvm(t,x,parms) Use:Function to calculate derivative of multispecies Lotka-Volterra equations Input: t:time (not used here,because there is no explicit time dependence) x:vector containing current abundance of all species parms:dummy variable,which is not used here (normally used to pass on parameter values,but not needed here because a and r are defined globally) Output: dx:derivative of Lotka-Volterra equations at point x 1vm <-function(t,x,parms){ dx<-(r-a%*%x)*x list(dx) Here,dr,r and x are vectors (whose length corresponds to the number of species)and a is an N x N matrix (where N is the number of species).The symbol %denotes matrix multiplication in R.(Note that x and da correspond to n and dn/dt in eq.3.)The command list at the end of the function just returns the values of dc. The numerical integration of the multi-species LVM can then be done using the function lsoda from the R library odesolve.A function to integrate a model such as the lvm looks for example like this: ### n.integrate(time,init.x,model) 3
aij > 1 then the effect of species j on i is stronger that the effect of species i on itself. If all aij > 0 then the multispecies LVM describes a community with only competitive interactions (i.e. without predator-prey relations between species). Predator-prey relations are characterized by interactions where aij 0. (Here species i is a predator of species j). Note, however, that for reasons of simplicity we will focus only on competitive interactions. (One of the advanced topics suggested below is to extend the approach to include predator-prey type interactions.) 2.1 Numerical integration of the multi-species LVM in R The accompanying starting R script provides an implementation for the numerical solution of the multi-species LVM. Please consult this script in parallel to reading this section. The R script is commented extensively. The numerical simulation of the multi-species LVM can be programmed in a few lines in R (if you know how!). Essentially what you need is a function that defines the derivative of the LVM (i.e. the right hand side of equation 3). This function looks like this: ### # lvm(t,x,parms) # Use: Function to calculate derivative of multispecies Lotka-Volterra equations # Input: # t: time (not used here, because there is no explicit time dependence) # x: vector containing current abundance of all species # parms: dummy variable, which is not used here (normally used to pass on # parameter values, but not needed here because a and r are defined globally) # Output: # dx: derivative of Lotka-Volterra equations at point x lvm <- function(t,x,parms){ dx <- (r - a %*% x) * x list(dx) } Here, dx, r and x are vectors (whose length corresponds to the number of species) and a is an N × N matrix (where N is the number of species). The symbol % ∗ % denotes matrix multiplication in R. (Note that x and dx correspond to n and dn/dt in eq. 3.) The command list at the end of the function just returns the values of dx. The numerical integration of the multi-species LVM can then be done using the function lsoda from the R library odesolve. A function to integrate a model such as the lvm looks for example like this: ### # n.integrate(time,init.x,model) 3
Use:Numerical integration of model Input: t:list with elements time$start,timeSend,and time$steps,giving start and endpoint of integration and the number of time points in between init.x:vector containing initial values (at time time$start)of all species model:name of the function to integrate (here lvm) Output: data frame with n+1 columns.The first column contains the time points at which x is evaluated.The next n columns are the values of the n species at these time points Description: Generates a vector of time points for the integration and uses function lsoda (from library odesolve)to integrate the model n.integrate <-function(time=time,init.x=init.x,model=model){ t.out <-seq(time$start,time$end,length=time$steps) as.data.frame(lsoda(init.x,t.out,model,parms=parms)) In the function we first generate a vector of times t.out between time$start and timeSend of length time$steps.Then we call the function lsoda,which requires initial values for the species vector x.These are specified by the vector init.x.The variable model stands for the model to integrate,which is here the function Ivm.The command as.data.frame returns the output of lsoda as an R dataframe. If we have specified all the relevant parameters and initial values then the LVM can be integrated by executing the following steps in R.First load the library odesolve by executing library(odesolve)on the R command line.This only needs to be done once at the beginning of an R session and it requires that the library odesolve is installed on your computer.Next you call n.integrate(time=time,init.x=init.x,model=1vm),where time is for example defined as time<-list(start=0,end=30,steps=100).Note,however,that you can only run the model provided you have specified all relevant parameters and initial values,such as the interaction matrix a,the growth rates r,the initial values for all species init.x,etc. 2.2 Generate and initialize parameter values The initialization of all parameters is in principle straightforward.For example you could execute the following commands in R: Number of Species n<-10 Generate n uniformly distributed random values (between 0 and 1)for the growth rate vector r <-runif(n) Generate n uniformly distributed random values for initial values of species
# Use: Numerical integration of model # Input: # t: list with elements time$start, time$end, and time$steps, giving start and # endpoint of integration and the number of time points in between # init.x: vector containing initial values (at time = time$start) of all species # model: name of the function to integrate (here lvm) # Output: # data frame with n+1 columns. The first column contains the time points at which # x is evaluated. The next n columns are the values of the n species at these # time points # Description: # Generates a vector of time points for the integration and uses function lsoda # (from library odesolve) to integrate the model n.integrate <- function(time=time,init.x= init.x,model=model){ t.out <- seq(time$start,time$end,length=time$steps) as.data.frame(lsoda(init.x,t.out,model,parms=parms)) } In the function we first generate a vector of times t.out between time$start and time$end of length time$steps. Then we call the function lsoda, which requires initial values for the species vector x. These are specified by the vector init.x. The variable model stands for the model to integrate, which is here the function lvm. The command as.data.frame returns the output of lsoda as an R dataframe. If we have specified all the relevant parameters and initial values then the LVM can be integrated by executing the following steps in R. First load the library odesolve by executing library(odesolve) on the R command line. This only needs to be done once at the beginning of an R session and it requires that the library odesolve is installed on your computer. Next you call n.integrate(time=time,init.x=init.x,model=lvm), where time is for example defined as time<- list(start=0,end=30,steps=100). Note, however, that you can only run the model provided you have specified all relevant parameters and initial values, such as the interaction matrix a, the growth rates r, the initial values for all species init.x, etc. 2.2 Generate and initialize parameter values The initialization of all parameters is in principle straightforward. For example you could execute the following commands in R: # Number of Species n<-10 # Generate n uniformly distributed random values (between 0 and 1) for the # growth rate vector r <- runif(n) # Generate n uniformly distributed random values for initial values of species 4
init.x <-runif(n) Generate n x n uniformly distributed random values for interaction matrix a <-matrix(runif(n*n),nrow=n) Integration window time<-list(start=0,end=30,steps=100) dummy variable for lvm()function defined above parms <-c(0)###dummy variable (can have any numerical value) Note that the dummy variable parms is required for the integration routine,but is of no further importance. Once you have specified these initial values you can numerically integrate the LVM by calling n.integrate(time=time,init.x=init.x,model=lvm).The output is a dataframe with the columns "time",for the time points,and "1","2",...for the abundances of the species at these time points.The start and end time and the time intervals in between are specified by the parameter vector time. 2.3 Notes on generating the interaction matrix There are several points that one needs to consider when generating the interaction matrix aij: The values ai have to be greater than zero,because otherwise a species might grow to infinity.Therefore,when generating the values for the interaction matrix one has to take care that the diagonal elements of the matrix (i.e.the ai)are greater than zero. When the interaction matrix aii is generated by drawing N x N uniformly distributed random numbers then essentially all species interact with all others (since the probability of generating a value of 0 vanishes).However,in natural ecosystems there will be many species that do not interact and therefore should have interaction parameters with value 0. 2.4 Tools to plot the output of the simulation The R script offers several functions to plot the generated output plot.lvm.time:This function uses the output of the function n.integrate to plot the time course of a simulation of the LVM. plot.matrix:This function plots the interaction matrix. plot.frequency:This function uses the output of the function n.integrate to plot the frequencies of all species at the end of the simulation. Examples for plots generated by these functions are given in figure 1.To understand how these functions work in detail,please refer to the R script file. 5
init.x <- runif(n) # Generate n x n uniformly distributed random values for interaction matrix a <- matrix(runif(n*n),nrow=n) # Integration window time<- list(start=0,end=30,steps=100) # dummy variable for lvm() function defined above parms <- c(0) ### dummy variable (can have any numerical value) Note that the dummy variable parms is required for the integration routine, but is of no further importance. Once you have specified these initial values you can numerically integrate the LVM by calling n.integrate(time=time,init.x=init.x,model=lvm). The output is a dataframe with the columns ”time”, for the time points, and ”1”, ”2”, . . . for the abundances of the species at these time points. The start and end time and the time intervals in between are specified by the parameter vector time. 2.3 Notes on generating the interaction matrix There are several points that one needs to consider when generating the interaction matrix aij : • The values aii have to be greater than zero, because otherwise a species might grow to infinity. Therefore, when generating the values for the interaction matrix one has to take care that the diagonal elements of the matrix (i.e. the aii) are greater than zero. • When the interaction matrix aij is generated by drawing N × N uniformly distributed random numbers then essentially all species interact with all others (since the probability of generating a value of 0 vanishes). However, in natural ecosystems there will be many species that do not interact and therefore should have interaction parameters with value 0. 2.4 Tools to plot the output of the simulation The R script offers several functions to plot the generated output. • plot.lvm.time: This function uses the output of the function n.integrate to plot the time course of a simulation of the LVM. • plot.matrix: This function plots the interaction matrix. • plot.frequency: This function uses the output of the function n.integrate to plot the frequencies of all species at the end of the simulation. Examples for plots generated by these functions are given in figure 1. To understand how these functions work in detail, please refer to the R script file. 5
8 10 1520 25 30 Sme 4 8 10 species index Figure 1:Plots generated by the plotting functions plot.lvm.time (top),plot.matrix(middle) and plot.frequency (bottom). 6
0 5 10 15 20 25 30 0.0 0.5 1.0 1.5 time abundance 1 2 3 4 5 6 7 8 9 10 10 9 8 7 6 5 4 3 2 1 2 4 6 8 10 0.0 0.5 1.0 1.5 species index abundance Figure 1: Plots generated by the plotting functions plot.lvm.time (top), plot.matrix (middle) and plot.frequency (bottom). 6
2.5 Repeated simulations The R script also offers the possibility to study the invasion of new species.In essence this is done by repeatedly running the LVM model,testing which species have gone extinct (i.e. have fallen in frequency below a pre-defined cutoff level),and replacing these extinct species by new ones that are allowed to invade into the ecosystem(see example loop in the R script).In order to replace the extinct species with new ones one needs to update the interaction matrix, the growth rate vector,and the init.x vector,for the initial values of the species at the start of the integration.The function cutoff can be used to determine which species went extinct and which survived at the end of a simulation.The function cutoff returns a vector with the indices of the extinct species.This vector is used by the function generate.parameters to update the growth rate vector and the interaction matrix.The function generate.parameters has another input parameter called sparse,which determines the fraction of all interactions that are nonzero.This parameter can take any value between 0 and 1.A value of sparse of 1 implies full connectivity.The function generate.parameters also makes sure that there are no 0 entries in the diagonal of the interaction matrix(i.e.all species inhibit their own growth). For details please refer to the R script. 3 Exercises 3.1 Basic exercises Ebl:How does ecosystem stability depend on size (i.e.the number of species present)?Hint: Start with a random N x N matrix,simulate the LVM and count how many species are present at the end of a (sufficiently long)simulation.Plot N against the number of surviving species. Eb2:How does ecosystem stability depend on connectivity (i.e.the fraction of non-zero entries in the interaction matrix)? Eb3:How can one define or measure ecosystem stability?Stability with regard to what?Think about measures that would be useful for natural systems.(This exercise is not primarily for simulation,but for discussion in order to think beyond the simulation. Eb4:Does the coexistence of a set of species depend on the order in which they are introduced into an ecosystem?Hint:Start with a given interaction matrix,but introduce species in different orders. 3.2 Advanced/additional exercises Eal:How does an ecosystem respond to the removal of a species?What is the average effect and what is the range of effects?Does the removal of an abundant species from the ecosystem have a stronger effect than the removal of a rare one?Does the removal of a species with high connectivity have a stronger effect than the removal of one with low connectivity?
2.5 Repeated simulations The R script also offers the possibility to study the invasion of new species. In essence this is done by repeatedly running the LVM model, testing which species have gone extinct (i.e. have fallen in frequency below a pre-defined cutoff level), and replacing these extinct species by new ones that are allowed to invade into the ecosystem (see example loop in the R script). In order to replace the extinct species with new ones one needs to update the interaction matrix, the growth rate vector, and the init.x vector, for the initial values of the species at the start of the integration. The function cutoff can be used to determine which species went extinct and which survived at the end of a simulation. The function cutoff returns a vector with the indices of the extinct species. This vector is used by the function generate.parameters to update the growth rate vector and the interaction matrix. The function generate.parameters has another input parameter called sparse, which determines the fraction of all interactions that are nonzero. This parameter can take any value between 0 and 1. A value of sparse of 1 implies full connectivity. The function generate.parameters also makes sure that there are no 0 entries in the diagonal of the interaction matrix (i.e. all species inhibit their own growth). For details please refer to the R script. 3 Exercises 3.1 Basic exercises Eb1: How does ecosystem stability depend on size (i.e. the number of species present)? Hint: Start with a random N × N matrix, simulate the LVM and count how many species are present at the end of a (sufficiently long) simulation. Plot N against the number of surviving species. Eb2: How does ecosystem stability depend on connectivity (i.e. the fraction of non-zero entries in the interaction matrix)? Eb3: How can one define or measure ecosystem stability? Stability with regard to what? Think about measures that would be useful for natural systems. (This exercise is not primarily for simulation, but for discussion in order to think beyond the simulation.) Eb4: Does the coexistence of a set of species depend on the order in which they are introduced into an ecosystem? Hint: Start with a given interaction matrix, but introduce species in different orders. 3.2 Advanced/additional exercises Ea1: How does an ecosystem respond to the removal of a species? What is the average effect and what is the range of effects? Does the removal of an abundant species from the ecosystem have a stronger effect than the removal of a rare one? Does the removal of a species with high connectivity have a stronger effect than the removal of one with low connectivity? 7
Ea2:How does an ecosystem respond to the invasion of a new species?Does invasion of a new species lead to the extinction of another species?How does the effect of invasion depend on the size of the ecosystem?How does it depend on the connectivity of the invading species and the connectivity among the other species? Ea3:Are "evolved"ecosystems more stable than random ones?Hint:Start with random interaction matrix and simulate the LVM until no further species goes extinct.Test the stability of the remaining system by randomly removing one of the species and determining how many additional species go extinct.Similarly,test "evolved"ecosystems for stability against extinction of a member species.By "evolved"system I mean a system that has been generated through multiple rounds of invasion of new species.Choose"evolved"and random systems of similar size. Ea4:How do the stability properties change if some interactions are predatory?This is com- putationally an advanced topic since it requires some new programming. 8
Ea2: How does an ecosystem respond to the invasion of a new species? Does invasion of a new species lead to the extinction of another species? How does the effect of invasion depend on the size of the ecosystem? How does it depend on the connectivity of the invading species and the connectivity among the other species? Ea3: Are ”evolved” ecosystems more stable than random ones? Hint: Start with random interaction matrix and simulate the LVM until no further species goes extinct. Test the stability of the remaining system by randomly removing one of the species and determining how many additional species go extinct. Similarly, test ”evolved” ecosystems for stability against extinction of a member species. By ”evolved” system I mean a system that has been generated through multiple rounds of invasion of new species. Choose ”evolved” and random systems of similar size. Ea4: How do the stability properties change if some interactions are predatory? This is computationally an advanced topic since it requires some new programming. 8