At this point I have been able to simulate two species predator-prey dynamics with and without density dependence. The next logical step is then to add another species, so that’s what I did.

library(deSolve) library(scatterplot3d) parameters<-c(r=2,alpha=.1,alpha2=.5,e=.1,e2=.1,q=.5,K=5000) state<-c(N=1000,N1=80,N2=20) lvmodel<-function(t,state,parameters){ with(as.list(c(state,parameters)),{ #rate of change dN<-r*N*(1-N/K)-alpha*N*N1 dN1<-e*alpha*N*N1-alpha2*N1*N2 dN2<-e2*alpha2*N1*N2-q*N2 #return rate of change list(c(dN,dN1,dN2)) }) } times<-seq(0,200,by=.1) out<-ode(y=state,times=times,func=lvmodel,parms=parameters) scatterplot3d(out[,2],out[,3],out[,4],xlab="N",ylab="N1",zlab="N2",typ="o",pch=20,cex.symbols=.75)

A quick note…

The equations for the dynamics of the basal and top species in this three species system are the same as those for the two species system. I was not completely sure what to do with the intermediate species. So, I assumed that the growth of the intermediate species was the numerical response to herbivory, while its mortality was a function of predation by the top species. That was what made sense to me, but there may be a better way to do it.

I visualized the output with the *scatterplot3d* package:

This was kind of cool, I was able to pick out some parameters that led to a stable equilibrium between the three species. Actually, one thing I noticed was that much as Hairston, Smith , and Slobodkin (HSS) argued, this world of three species I created is green.

For those unfamiliar with the “world is green” hypothesis, a paper by HSS made the argument that there tends to be an abundance of plants as a result of top down control. basically, predators keep herbivore populations small, allowing plants to grow without being eaten. The paper is very well written, but in my opinion is a little heavy in argument and light on data.

Regardless, back to my simulated system. The reason I bring up HSS is that I found that the equilibrium abundances of my three species, given the parameters I selected, matched HSS predictions (abundant plants, small herbivore population, larger predator population). Of course that doesn’t necessarily mean anything, I just thought it was a fun coincidence.

I played a little bit with the parameters of the above model, and got a few different trajectories.

I think this is one way to model three species Lotka-Volterra systems, but doing this would be a remarkably inefficient way to scale up to more than three, or even anything more complex than A eats B eats C. So that is my next move, multispecies Lotka-Volterra equations in R.

One thing that I think will help me to figure out how to do this is a section from a class offered at the Swiss Federal Institute of Technology Zurich on modelling in Population and Evolutionary Biology. They have a section on Stability and Complexity of Model Ecosystems. Along with the reading they have provided some sample code:

</p> <p style="text-align: left;">###################################################################### ### Stability of ecosystems: Lotka-Volterra Multispecies model ###################################################################### ################################### # FUNCTION DEFINITIONS ################################### ### # 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 the 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) } ### # n.integrate(time,init.x,model) # 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)) } ### # cutoff(out,tol) # Use: Get index of species (i) that are extinct (i.e. whose frequency is smaller than tol) # (ii) species that survive (i.e. whose freq is larger than tol) and (iii) return a vector # of the abundance of species, in which those with freq < tol are set to 0. # Input: # out: This is the output of n.integrate (i.e. a data frame with a column for the # time points and the abundance of the n species) # tol: The threshold value below which a species is considered to be extinct # Output: # list with three entries # (i) indices of species that are extinct # (ii) indices of species that survive # (iii) vector of species abundances where those with freq < tol are set to 0. (This vector # can be used as the new initial conditions to continue the integration of the model # Description: # Takes output of n.integrate and determines which species # are more abundant than tol at the last time point. cutoff <- function(out,tol){ n.end <- length(out[,1]) frequency <- out[n.end,2:(n+1)] #Note that species abundance is in columns 2 to n+1 extinct <- which(frequency<tol) survive <- which(frequency>=tol) frequency[extinct] <- 0 list(extinct=extinct,survive=survive,freq=as.numeric(frequency)) } ### # generate.parameters(mat,rep,index,sparse) # Use: Gets reproduction rate and interaction matrix and generates new reproduction rates and interaction # coefficients of species with given set of indices # Input: # mat: current matrix of interaction coefficients # rep: current reproduction rates of the species # index: vector of indices of species for which new parameters have to be created (use for # example output of cutoff() as indices for extinct species) # sparse: Sparsity of interaction matrix. The parameters sparse controls how many # interactions are set to 0. sparse has to between 0 and 1, where sparse = 1 # corresponds to full connectivity. # Output: # list containing updated reproduction rates and interaction matrix generate.parameters <- function(mat,rep,index,sparse){ length.index <- length(index) # generate random reproduction rates for species with indices in index vector rep[index]=runif(length.index) # generate uniformly distributed random interaction coefficients to replace rows of # interaction matrix with corresponding indices m1 <- matrix(runif(n*length.index),ncol=n) # generate matrix of 0's and 1's as a function of sparse (using rbinom()) m2 <- matrix(rbinom(n*length.index,1,sparse),ncol=n) # replace the corresponding rows in the interaction matrix mat[index,] <- m1*m2 # Do the same for the columns mat[,index] <- matrix(runif(n*length.index),nrow=n)*matrix(rbinom(n*length.index,1,sparse),nrow=n) # make sure that diagonal elements are not zero (i.e. regenerate random coefficients for diag elements) diag(mat)[index] <- runif(length.index) list(matrix=mat,reproduction=rep) } ### Plotting routines ### # plot.lvm.time(out) # Use: Plots all species against time # Input: # out: output of n.integrate plot.lvm.time <- function(out){ n <- length(out[1,])-1 t.range <- range(out[,1]) y.range <- range(out[,2:(n+1)]) plot(t.range,y.range,xlab="time",ylab="abundance",type="n") for (i in c(2:(n+1))) {lines(out[,1],out[,i],col=i)} } ### # plot.matrix(matrix,index) # Use: Plot interaction matrix (setting all species in index to zero) # Input: # matrix: interaction matrix # index: indices of species that are set to 0. (These would be the species # indices determined by function cutoff() that are set to zero because # the species abundance is smaller than tol). plot.matrix <- function(matrix,index){ m <- matrix m[index,] <-0 m[,index] <-0 image(t(m),axes=F,zlim=c(0,1),col=gray(c(32:0)/32),xlab="species affected",ylab="species affecting") axis(1, at = seq(0,1,length=n),labels = 1:n) axis(2, at = seq(0,1,length=n),labels = 1:n) } ### # plot.frequency(out) # Use: Plots abundance of species at the end of simulation # Input: # out: output of n.integrate() # Plot frequency at end of simulation plot.frequency <- function(out) { plot(c(1:n), out[length(out[,1]),2:(n+1)],type="h",xlab="species index",ylab="abundance") } ########################### # MAIN PROGRAM ########################### ### LOAD LIBRARIES #load R library for ordinary differential equation solvers library(deSolve) ### INITIALIZE PARAMETER SETTINGS # 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) # Number of Species n<-10 # Sparsity of interaction matrix (0<sparse<1) [sparse = 1 ==> full connectivity] sparse <- 0.2 # Generate parameters # first generate reproduction vector and interaction matrix containing only zeros. r <- numeric(n) a <- matrix(0,nrow=n,ncol=n) # Now generate random numbers for parameters according to function generate.parameters tmp <- generate.parameters(a,r,c(1:n),sparse) r <- tmp$reproduction a <- tmp$matrix # Initial conditions (i.e. frequencies of n species at time = time$start) init.x <- runif(n) # Threshold at which species are assumed to be extinct tol <- 0.0001 ### Example for a loop that integrates the LVM and allows new species to invade outt <- init.x for (i in c(1:20)){ # integrate lvm model out <-n.integrate(time=time,init.x=init.x,model=lvm) # plot time course plot.lvm.time(out) # save out and attach it to outt outt <- rbind(outt,out) # determine which species went extinct and which survived tmp <- cutoff(out,tol) # save index of extinct species (i.e. species with frequency smaller than tol) index.extinct <- tmp$extinct # generate new parameters (reproduction rates and interaction coefficients) to replace extinct species tmp2 <- generate.parameters(a,r,index.extinct,sparse) a <- tmp2$matrix r <- tmp2$reproduction # generate new init.x to continue integration init.x <- tmp$freq # introduce new species with a frequency of 10* tol init.x[index.extinct] <- 10* tol # generate new time window to continue integration time <- list(start=time$end, end = time$end+time$end-time$start, steps=100) } # plot entire time course plot.lvm.time(outt) </p> <p style="text-align: left;">

I found this today, and it is pretty cool. I plan on going through this with a fine tooth comb to understand all the different functions and how they work together. I hope to adapt this approach to help me with my own research.

OnOne

–>