Two species predator prey systems with R

It has taken me entirely too long to learn how to use R to solve differential equations. Clearly, with my interest in modelling dynamics on food webs, being able to solve differential equations is critical. I had known about the package odesolve, but apparently that has been replaced with deSolve. Turns out, its fairly easy to model a two species predator prey system. Unfortunately, I am still figuring out how to work out parameterization. My Lotka-Volterra predator prey equations give what I imagine is a strange output (which may most likely be the result of my using strange values for the parameters), but when I include prey density dependent growth it seems to be more reasonable.

Screen shot 2013-01-31 at 12.11.01 PM

This is the code for modeling a two species Lotka-Volterra system that I am using. The values for the parameters are completely fabricated (not based on anything other than whimsy). I do wonder whether the code above can reproduce published data, but I have to find some published data first (although I suppose that really shouldn’t be too hard). When I run the above code (also posted in copy/paste-able form below) I get the following result:

This is the graphical output of the code that is posted above.

This is the graphical output of the code that is posted above.

My understanding is that this system of equations should give me stable limit cycles around some fixed equilibrium point. Near as I can tell this is close, but seems to be a little wonky. I am thinking it could be unrealistic parameters but my expectation was that I would get something a little more circular. I have, with other values, gotten a result to be somewhat egg shaped which I thought was pretty cool.

Nonetheless I am going to keep playing around with this basic code form and using deSolve. My goal is to move from this two species system to increasing numbers of species (3, 4, 5, 6) as part of a research project I am working on for understanding stability of trophic chains of varying length (more on that later). From there I want to move on to n-species systems with changing parameter values.

…Updated since I just learned how to make the r code more presentable, thanks to Chit Chat R


library(deSolve)

parameters<-c(r=1.15,alpha=.01,q=.4)
 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<-(alpha*N*P)-(q*P)

list(c(dN,dP))
 })
 }

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=.4/.01)
 abline(h=1.15/.01)

Advertisements
This entry was posted in Research and tagged , , , , , . Bookmark the permalink.

9 Responses to Two species predator prey systems with R

  1. Steve says:

    Hi there, looks good so far, the lotka-volterra equations are a great model system to play with! You note that your output looks a bit strange. By decreasing the stepsize to something like:
    times<-seq(0,100,by=0.1)
    should make your phase space diagram appear a lot less wonky!

  2. Mike Fowler says:

    Hi Jon,
    Another common trick is simply to present the dots without connecting lines, and/or present the data as log(Abundance). This might give a more ‘circular’ cycles.

    This is a great post – I also keep meaning to invest more effort into switching to R and this post give a useful ‘crib’ sheet!

    I’d be interested to know how long your simulation code takes to run, e.g., for a comparison with MatLab. Also, what sort of integration does the R solver use? Is it e.g., 2nd-3rd order, 4th-5th, or something else? You seem to force the output in time steps of 1 unit (or 0.1 as Steve suggests above), but this might not be the step size the solver is actually using for integration.

    • Jon Borrelli says:

      I admit, I had thought about just eliminating the lines, but I feel it is more aesthetically pleasing with them. I will have to see what happens when I use log abundance, that is an interesting idea.
      In terms of the system time:
      > system.time(
      + out
      That is how long it took to run the code from the post. I have also tried it with up to 10000 time steps and it takes just a little longer (although I suppose it could really ad up if you are iterating the model multiple times).

      I am not 100% certain how exactly the method of integration works. I am still working on learning that. Although, from what I have read the ode function in deSolve allows the user to specify the method of integration from a relatively large list. The default is lsoda, but there are also a number of Runge-Kutta methods (if that means anything to you). There are a few nice documents out there about the different methods of integration used in this package, found most easily by googling “solving differential equations in r.”

      Lastly, that is an interesting point about the time steps. To my understanding the solver can either be calculating based on the time steps provided, or based on an independent set of time steps and then interpolating the results. I think that which method is used is dependent on which method of integration is specified.

      I am going to have to do a bit more reading on this I think. Thanks for posing such interesting questions!

      • Jon Borrelli says:

        Just noticed it didn’t paste the actual system time…

        > system.time(
        + out<-ode(y=state,times=times,func=lvmodel,parms=parameters)
        + )
        user system elapsed
        0.222 0.002 0.223

  3. Ben Bolker says:

    (1) the graphs are also prettier (I think) on a log scale — your starting conditions are such that the limit cycle gets down to prey densities of around 10^{-4}. (2) most of the deSolve integrators have adaptive step sizes (see vignette(“deSolve”), section 2).

    • Jon Borrelli says:

      Thanks! I tried out using log scale and it works very nicely. That is something I really should have thought of. That’s a good tip about the time step, thanks for pointing out that section. Now things are starting to get a little clearer for me.

  4. Jeremy Fox says:

    If you’re interested, Hank Stevens has a good book on teaching theoretical population ecology with R. It’s got R code for things like standard continuous time predator prey models, competition models, and much, much more.

    http://www.springer.com/life+sciences/ecology/book/978-0-387-89881-0

    • Jon Borrelli says:

      Thanks for pointing that out Jeremy. I had been looking for Use R! books, but I had not seen this one. I am going to have to take a look into getting this book, it would be great to have a book that could teach me the basics of using these models in R without my having to stumble through it on my own!

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s