Stage 1:

Recently came across this concept while talking to my roommate. A first glance through the Wikipedia page showed that it was basically a mathematical model called the Lotka Volterra Equation . The Programmer in me decided to implement it in R, and the cynic in me told me to Google it first. The first search gave the following results:

Now, I was playing around with this code, which is as follows,

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)

I must say that this code works perfectly fine. The only problem was I have to manually change the \(\alpha\) , \(\beta\) , \(\gamma\) and \(\delta\) values manually to see the plots.

Idea : Why not make a function out of it?

I made some minor changes here and there to get the following code in the form of a function:

PrPred <- function(a,b,g,d){
  library(deSolve)

  Pars <- c(a, b, g, d)
  State <- c(x = 10, y = 10)


  LotVmod <- function (Time, State, Pars) {
    with(as.list(c(State, Pars)), {
      dx = x*(a - b*y)
      dy = -y*(g - d*x)
      return(list(c(dx, dy)))
    })
  }

  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("Prey", "Predator"), lty = c(1,2), col = c(1,2), box.lwd = 0)

}

In the above code \(a = \alpha\), \(b= \beta\), \(g = \gamma\) and \(d= \delta\).

Now, I remembered something I read on Hilary Parker’s blog|Make an R package from scratch, where she says,

Hadley Wickham gets full credit for envisioning that R packages should be the easiest way to share code, and making functions/resources that make it so easy to do so.

So, I decided to make a R Package out of it, and I did. Find it on Github|Prey-Predator-Model.

Stage 2:

When I was changing the \(\alpha\) , \(\beta\) , \(\gamma\) and \(\delta\) values and playing around the plots, I could not understand much of it. There was some rigorous mathematics going on in this concept.

The \(\alpha\) , \(\beta\) , \(\gamma\) and \(\delta\) have distinct meanings. They are as follows:

\(\alpha\) - The growth rate of Prey,
\(\beta\) - Rate at which Predators destroy Prey,
\(\gamma\) - The death rate of predators,
\(\delta\) - The rate at which predators increase by consuming prey.

For more refer : Lotka-Volterra Equations – On Wolfram Mathworld.

The equations that govern the Prey-Predator model are as follows:

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

Here,

\(x\) = Number of Prey.
\(y\) = Number of Predator

Now, we try to break equations into pieces and analyse each term separately,

Equation 1:

A prey population (\(x\)) increases at a rate \(dx = \alpha xdt\) (proportional to the number of prey).
A prey population is simultaneously killed at a rate \(dx = -\beta x y dt\)(proportional to the product of the numbers of prey and predators).

Equation 2

A predator population (\(y\)) decreases at a rate \(dy=-\gamma ydt\) (proportional to the number of predators).
And increases at a rate \(dy=\delta xydt\) (proportional to the product of the numbers of prey and predators).

N.B. There are multiple values for which the model doesn’t hold good, the function that I have written will help you to see the behavior of the Population of Prey and Predator with time as you alter \(\alpha\), \(\beta\), \(\gamma\) and \(\delta\).

References :

The Prey Predator model

Stage 1:

Stage 2: