Lotka-Volterra model ~ intro

Posted on March 30, 2010 by apeescape in R bloggers | 0 Comments

[This article was first published on mind of a Markov chain » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

So many know about the Lotka-Volterra model (i.e. the predator-prey model) in ecology. This model portrays two species, the predator (y) and the prey (x), interacting each other in limited space.

$frac{dx}{dt} = x(alpha - beta y)\ frac{dy}{dt}=-y(gamma - delta x)$

The prey grows at a linear rate ( $alpha$ ) and gets eaten by the predator at the rate of ( $beta$ ). The predator gains a certain amount vitality by eating the prey at a rate ( $delta$ ), while dying off at another rate ( $gamma$ ).

Given this base, we can ask questions like, what parameterizations can we expect to find a coexistence between the fox and the hare (for example)?

Let’s choose some values for the model: $alpha = 2, beta = .5, gamma = .2, delta = .6$ . These values assume a weaker growth of the rabbits relative to the strength of the death of foxes. Below, I simulated these values in R.

And we get coexistence, they live happily forever after. With this simple model, we can play around by generalizing (logistic growth of prey, etc.). I will put up some posts doing so.

The way to do this in R is as follows (just use the deSolve package, which will supersede the odesolve package):

 library(deSolve)

LotVmod <- function (Time, State, Pars) {
	with(as.list(c(State, Pars)), {
		dx = x*(alpha - beta*y)
		dy = -y*(gamma - delta*x)
		return(list(c(dx, dy)))
	})
}

Pars <- c(alpha = 2, beta = .5, gamma = .2, delta = .6)
State <- c(x = 10, y = 10)
Time <- seq(0, 100, by = 1)

out <- as.data.frame(ode(func = LotVmod, y = State, parms = Pars, times = Time))

matplot(out[,-1], type = "l", xlab = "time", ylab = "population")
legend("topright", c("Cute bunnies", "Rabid foxes"), lty = c(1,2), col = c(1,2), box.lwd = 0)

Filed under: deSolve, Food Web, ODEs, R