Correct use of deSolve in ecological modelling of a predator prey system

Question

I have a predator-prey model with parameters and initial values as specified. I solve the differential equations two ways here 1. using a for loop and 2. using the deSolve package. I believe the for loop is correct, and should give output as seen in the plot below.

##################
#For loop attempt
##################
r=0.1; k=100; v=40; s=.1; tbeg=0; tend=1200; nints=1200
c=.06; a=.12; predator0=c(15); prey0=c(50)

dt=(tend-tbeg)/nints
prey=matrix(0,nints+1,length(prey0))
predator=matrix(0, nints+1, length(predator0))
predator[1, ]=predator0
time=numeric(nints+1)
prey[1, ]=prey0
for(i in 1:nints) {
    dprey=r*prey[i, ]-(r*prey[i, ]*prey[i, ]/k)-(s*prey[i, ]*predator[i, ])/(v+prey[i, ])*dt
    prey[i+1, ]=prey[i, ]+dprey
    dpredator=(a*prey[i, ]*predator[i, ])/(v+prey[i, ])-(c*predator[i, ])*dt
    predator[i+1, ]=predator[i, ]+dpredator
    time[i+1]=time[i]+dt}
matplot(time, prey, type="l", lty=1, main="Case 1:  Predator and Prey Populations over Time", ylab="Population", xlab="Time")
points(time, predator, type="l", lty=2)

I am attempting to use the package deSove to solve this system of DE's and expect similar output. The code runs, but provides a different answer. (Following "A Primer on Ecology" examples)

##################
#deSolve attempt 
##################
library(deSolve)
#the case of a prey with a logistic growth and predator functional response
predprey_FuncResp  <- function(t, y, parms) {
 n0 <- y[[1]]
 p0 <- y[[2]]
 with(as.list(parms), {
    dpdt <- (a*n0*p0)/(v+n0) - c*p0
#    dndt <- r*n0 - (r*n0^2)/k - (s*n0*p0)/(v+n0)
    dndt <- r*n0*(1-n0/k) - (s*n0*p0)/(v+n0)
    return(list(c(dpdt, dndt)))
    })
 }
parms <- c(a=.12, c=.06 , r=.1,s=.1,k=100, v=40)
 Tmax = 1000 # time horizon  
 TimeStep = 1 # integration time step
 Time <- seq(0, Tmax, by = TimeStep)  # the corresponding vector
LV.out <- lsoda(c(n0 = 50, p0 = 15), Time, predprey_FuncResp, parms)
matplot(LV.out[,1],LV.out[,-1], type='l')

I assume I am using deSolve incorrectly but cannot see my error. Thanks to anyone who takes the time to look at this.

Answer 1

Your problem is that you have the order of the terms switched when you return them from the gradient function -- you need list(c(dndt, dpdt)) if your initial conditions are {n0,p0}.

I made a few additional cosmetic tweaks.

library(deSolve)
predprey_FuncResp  <- function(t, y, parms) {
    with(as.list(c(y,parms)), {
         dpdt <- (a*n0*p0)/(v+n0) - c*p0
         dndt <- r*n0*(1-n0/k) - (s*n0*p0)/(v+n0)
         return(list(c(dndt, dpdt)))
    })
}
parms <- c(a=.12, c=.06 , r=.1,s=.1,k=100, v=40)
Tmax = 1200 # time horizon  
TimeStep = 1 # integration time step
Time <- seq(0, Tmax, by = TimeStep)
LV.out <- ode(c(n0 = 50, p0 = 15), Time, predprey_FuncResp,parms)
par(las=1,bty="l")
matplot(LV.out[,1],LV.out[,-1], type='l', xlab="time", ylab="density")