Examples of ODE's in R from classic ecological models with a simple M-W extension on I

ecoODE.R

#' Description: Code using deSolve for some simple ecological ODE models
#' There are 3 examples, one is a simple LV Pred-Prey example with plots
#' The next is a Macarthur-Wilson model with plots of both the species curve and the equlibrium points
#' The last plot is a quick crack at putting stochasticity into the model by adding oscillations in I with t.
#' LV code is modified from a blog post here: http://assemblingnetwork.wordpress.com/2013/01/31/two-species-predator-prey-systems-with-r/
#' Other additions are my own
#' 
#' Author: Edmund Hart
#' Date: 3/28/2013
#' E-mail: [email protected]



library(deSolve)

### Example using LV- Pred-Prey model
state<-c(N=500,P=20)

lvmodel<-function(t,state,parameters){
  with(as.list(c(state,parameters)),{
    #rate of change
    dN<-(r*N)-(alpha*N*P)
    dP<-(beta*N*P)-(q*P)
    
    list(c(dN,dP))
  })
}
parameters<-c(r=1.5,alpha=.01,beta=.004,q=.8)

times<-seq(0,100,by=1)
out<-ode(y=state,times=times,func=lvmodel,parms=parameters)
plot(out[,2],out[,3],typ="o",ylab="Predator Abundance",xlab="Prey Abundance")
abline(v=parameters[4]/parameters[3],lty=2,lwd=2)
abline(h=parameters[1]/parameters[2],lwd=2)

plot(out[,2],type='l', main="Predator-prey cycles from LV",xlab="Time",ylab="N")
lines(out[,3],col=2)




### Working with Macarthur-Wilson

state<-c(S=10)

mwmodel<-function(t,state,parameters){
  with(as.list(c(state,parameters)),{
    #Change in Species, follows the formulation in Gotelli 2001
    lambda <- I -(I/P)*S
    mu <- (E/P)*S 
    dS<- lambda - mu
      
    list(c(dS))
  })
}
parameters<-c(P=100,I=5,E=5)

times<-seq(0,100,by=1)
out<-ode(y=state,times=times,func=mwmodel,parms=parameters)
plot(out,xlab="Time",ylab="Number of Species", main="McArthur-Wilson model")

#### For plotting purposes
lambda <- function(I,P,S){return(I-(I/P)*S)}
mu <- function(E,P,S){return((E/P)*S)}
S <- 0:100
plot(S,lambda(parameters[2],parameters[1],S),type='l',ylab="Immigration rate",xlab="Number of species",main="McArthur-Wilson")
lines(S,mu(parameters[3],parameters[1],S),col=2)




####Here's a modification of MW to include cycling in imigration rates with time.
### Increasing the parameter M in the model will increase the amplitude.
### set to 1 for small effects.
state<-c(S=10,I=5)

modmwmodel<-function(t,state,parameters){
  with(as.list(c(state,parameters)),{
    #Change in Species, follows the formulation in Gotelli 2001
     
    mu <- (E/P)*S 
    dS<- I -(I/P)*S - mu
    dI <- cos(t)*M
    
    list(c(dS,dI))
  })
}
parameters<-c(P=100,E=5,M=10)

times<-seq(0,100,by=1)
out<-ode(y=state,times=times,func=modmwmodel,parms=parameters)
plot(out[,1],out[,2],xlab="Time",ylab="Number of Species", main="McArthur-Wilson model",type='l')