This page will teach you how to build, in R, a model simulating a consumer-resource interaction.

We will create a simple difference equation population interaction model. We will pretend it represents a population of parasitoids which develop in and kill individuals in a population of hosts.

You will learn how to use a loop function (the for() function in R), as well as the basics of difference equation model building. I encourage you to play around with the parameters presented below to see how they change model output. Importantly, regarding the R code, there is more than one way to skin a cat. You may find a more elegant or computationally less intensive way of doing things. What I present below is simply what I find to be the most intuitive and helpful from a heuristic standpoint.


Let’s get started…

Before we get into a more complicated interaction model, let’s start by modelling a single host population with exponential growth.

Let’s say we want the population to grow so that \[H_{t+1} = H_te^r\] where, \(H\) is the host population, \(r\) is the intrinsic rate of increase, \(t\) is the time step, and \(e\) is a constant approximately equal to 2.7.

What this model is saying is that at the time point \(t+1\), the population density will be \(e^r\) times bigger than it was at the previous time point \(t\).

(Note: OK, I said we would model exponential growth, but technically this should be called “geometric growth” because it is a discrete time model. That is, time is modeled as fixed time points rather than continuously.)

To run a simulation of this model in R, we will create 5 things.

1. We will create an empty vector for the the host population, \(H\). We want this so we can later fill in the elements for the NULL vector. Each element will represent the population density for each time step in the model.

H=NULL

2. If this population is to grow, we need a growth rate.

r=0.25

This parameter, \(r\), is called the intrisic rate of increase for the population. At each time step the population will grow at a rate of \(e^r\).

If you’re curious why we use the constant \(e\) instead of some other number, I think this guy explains it really well.

Sometimes you will see the above model written as
\[H_{t+1} = \lambda H_t\]

In this case, \(\lambda = e^r\). Conversely, \(r = log_e \lambda\). Because of the manner in which I will add subsequent pieces to this model, below, I find it easier to write the model using \(e^r\), instead of \(\lambda\).

3. We also need to tell R what the initial conditions are. In other words, what does the population density look like at time \(t=1\)? Let’s make a starting density for \(H\).

H[1]=5  

Now the first element of the vector H is 5. Simple.

4. We need to decide for how long we should run our model simulation. ….How about 50 time steps?

time=50

Piece of cake.

Now for the tricky part.

5. We’ve defined all of our parameters, our initial conditions, and how long we want to run the simulation for. Now we need to code the model itself. To do this, we need to find a way to make R calculate a future time step based on what happened in the previous time step. Remember our model in mathematical terms: \[H_{t+1} = H_te^r\] We want things that happen at time \(t+1\) to be dependent on what happened at time \(t\).

To do this, we will use a for() loop function. We’ll learn how this works below.

First, take a look at what I’ve got here.

for(i in 1:time){
  H[i+1]=H[i]*exp(r)
}

That right there is a geometric growth model.

This is how it works. The for() loop will run through the vector 1:time. Remember we set time=50 above? So 1:time is a vector that is all the integers in sequence from 1 to 50. It will represent \(t=1\), \(t=2\), \(t=3\), …, \(t=50\).

The for() loop, will interatively go through the operations in the {} brackets for as many elements as there are in the vector 1:time (i.e. 50). On the \(i^{th}\) interation, i will be equal to the \(i^{th}\) element in the vector 1:time. For example, in the first iteration, we have i=1, so R will compute \[H_{1+1} = H_1e^r\]

In R, this would be like writing H[1+1]=H[1]*exp(r). The command exp() raises \(e\) to whatever we put in the parenthesis, so exp(r) equals \(e^r\).

For the second iteration, i=2, so R will compute \[H_{2+1} = H_2e^r\]

The vector H, for which we had previously only defined its first element (i.e.H[1]=5), will get elements filled in each time the loop runs through an iteration.

Do you see how the for() command works? Really try to understand this part. for() loops are super useful in running simulations, doing bootstap analyses, and all sorts of other things in R. This is worth your time to learn.

Okay

Let’s take a look at all of our R code we’ve created thus far.

H=NULL

r=0.25

time=50

H[1]=5

for(i in 1:time){
  H[i+1]=H[i]*exp(r)
}

After we run this code, let’s take a look what H looks like.

We’ll plot H on the y-axis and time (or generations) on the x-axis. I won’t expain plot() commands here. I think the help menu (to view type ?plot)is suffient for the level of graphical display we’ll use here.

plot(1:(time+1),H,xlab="generations",ylab="host density")
lines(1:(time+1),H)

Sweet - geometric growth.

But how often do we see limitless exponential growth in an insect population? Let’s put an upper bound (or a carrying capacity) on that growth and make it logistic growth.

\[H_{t+1} = H_te^{r(1-H_t/K)}\]

What we’ve added here is a carrying capacity, \(K\). As the value of \(H\) approaches the value of \(K\), the effect of the intrinsic growth rate, \(r\), is modified by the term \((1-H/K)\).

If \(H\) is very low, this term will roughly equal \(1\), because \(1-tiny\) is approximately 1, so the growth rate will approximate the maximal growth rate (i.e. \(e^r\)).

If \(H\) is equal to \(K\), population growth will stop because \(1-1=0\), and \(e^0=1\).

If \(H\) overshoots \(K\), the term will become negative, because \(1-large=negative\), and \(e^{negative}\) is less than 1. This will cause the population to decline.

Make sense?

Let’s code the logistic growth model in R.

H=NULL #recreate H as NULL
H[1]=5 #define starting conditions

time=50 #how long are we running the simulation
r=0.25 #intrinsic growth rate
K=200 # carrying capacity

for(i in 1:time){
  H[i+1]=H[i]*exp(r*(1-H[i]/K))
}

Now that we have added carrying capacity to our R code, let’s see what our population looks like through time.

plot(1:(time+1),H,xlab="generations",ylab="host density")
lines(1:(time+1),H)

Cool.

Now the host population doesn’t just ridiculously grow to infinity as it did in the geometric growth equation above. With the carrying capacity parameter added, the population growth rate has become density-dependent. That means that the population growth rate depends on the density of the population. Specifically, the growth rate at \(t+1\) depends on the population density at \(t\).

In our earlier geometric growth model, while the population density depended numerically on the density of the previous time step, the growth rate itself was constant.

One thing to notice about difference equation models is that there is a lag in how the dynamics of the model responds to density. The lag occurs because the parameters are adjusted based on events that happened at a previous time step. This means there is a potential for the population to overshoot boundaries like the carrying capacity. For example, look what happens when we make the intrinsic growth rate larger.

H=NULL #recreate H as NULL
H[1]=5 #define starting conditions

time=50 #how long are we running the simulation
r=3 # larger intrinsic growth rate
K=200 # carrying capacity

for(i in 1:time){
  H[i+1]=H[i]*exp(r*(1-H[i]/K))
}

plot(1:(time+1),H,xlab="generations",ylab="host density")
lines(1:(time+1),H)

The host population overshoots its carrying capacity and oscillates back and forth around this value. I encourage you to play around with the parameters \(K\) and \(r\), as well as the initial host density, so that you get a better feel for how the model behaves.


Now let’s bring in some parasitoids.

The basic idea in a consumer-resource interaction model is to make one or both of the populations dependent on some aspect of the other population. In our example below, we will make growth of the parasitoid population dependent on the the density of hosts, and growth of the host population will be affected by how many parasitoids are around to attack them.

We’ll roughly follow Nicholson and Bailey and make parasitoid growth dependent on the number of hosts in the environment that it is able to attack.

\[H_{t+1} = H_te^{r(1-H_t/K)-aP_t}\] \[P_{t+1} = H_t (1-e^{-aP_t})\]

Okay, so what have we done here?

For starters we added an equation for \(P\) which explains the growth of the parasitoid population. Similar to the host population, the density at time \(t+1\) depends on what happened at time \(t\).

Another thing we did was add an important new parameter, a. This parameter is known variously as search efficiency or search area or sometimes just attack. Don’t call this parameter attack rate because it can’t actually be interepreted as a rate of attack (i.e. the number of hosts attacked per unit time), although it does directly affect the rate of attack. That attack rate for this equation can easily be calculated for a given time step, but that is beyond the scope of this lesson.

In the parasitoid equation, the term \((1-e^{-aP_t})\) is a proportion (0 to 1), and it tells us the proportion of hosts that are attacked. In the parasitoid equation, this value is multiplied by \(H_t\), which means some proportion of the host population is attacked, and those attacked hosts develop into parasitoids at time \(t+1\).

Notice that in the host equation we added to the \(e\) exponent the term \(-aP_t\), while in the parasitoid equation we have the term \((1-e^{-aP_t})\). If \((1-e^{-aP_t})\) is the proportion of hosts that get attacked, then \(e^{-aP_t}\) must be the proportion of hosts that don’t get attacked. Thus, \(e^{-aP_t}\) is often referred to as an escape term, or if it’s written as its own function, we call it the escape function.

So, in the host population, the growth rate at a given time point is determined by three things: 1) the proportion of hosts that escape parasitism, \(e^{-aP_t}\); 2) the carrying capacity rate-adjustment term \({(1-H_t/K)}\); and 3) the maximal growth rate, \(e^r\).

Make sense?

Okay, but what does this actually look like?

Here’s the R code and a plot of the dynamics through time.

P=NULL  #create an empty set for the parasitoids
H=NULL  #create an empty set for the hosts
P[1]=1  #define initial conditions for the parasitoids
H[1]=5  #define initial conditions for the hosts

time=50  #how many time steps to run our simulation
r=0.25  #intrinsic growth rate of the host
a=0.25  #search efficiency of the parasitoid
K=200  #carrying capacity of the host

#here's our model
for(i in 1:time){
  H[i+1]=H[i]*exp(r*(1-H[i]/K)-a*P[i])
  P[i+1]=H[i]*(1-exp(-a*P[i]))
}

#this code below will plot P and H on the y-axis, and time on the x-axis
plot(1:(time+1),H,xlab="generations",ylab="density",ylim=c(0,max(H)))
lines(1:(time+1),H,pch=1)
lines(1:(time+1),P,col="blue")
points(1:(time+1),P,pch=2,col="blue")
legend(.8*time,.9*max(c(H,P)),legend=c("H","P"),col=c("black","blue"),lty=c(1,1),pch=c(1,2))

That looks so sciencey.

The coding is a fairly straightforward addition to our previous logisitc model.

But what is our simulation telling us about the model???

It looks like the parasitoid and host populations are oscillating back and forth. But I wonder what happens further out. I don’t think 50 time steps is sufficient for us to figure out what’s going to happen here.

Any guesses what would happen if run the simulation longer?

Let’s find out. We’ll run it for 150 generations.

P=NULL  #create an empty set for the parasitoids
H=NULL  #create an empty set for the hosts
P[1]=1  #define initial conditions for the parasitoids
H[1]=5  #define initial conditions for the hosts

time=150  #how many time steps to run our simulation
r=0.25  #intrinsic growth rate of the host
a=0.25  #search efficiency of the parasitoid
K=200  #carrying capacity of the host

#here's our model
for(i in 1:time){
  H[i+1]=H[i]*exp(r*(1-H[i]/K)-a*P[i])
  P[i+1]=H[i]*(1-exp(-a*P[i]))
}

#this code below will plot P and H on the y-axis, and time on the x-axis
plot(1:(time+1),H,xlab="generations",ylab="density",ylim=c(0,max(H)))
lines(1:(time+1),H,pch=1)
lines(1:(time+1),P,col="blue")
points(1:(time+1),P,pch=2,col="blue")
legend(.8*time,.9*max(c(H,P)),legend=c("H","P"),col=c("black","blue"),lty=c(1,1),pch=c(1,2))

Huh?!

It looks like our populations oscillated back and forth, spinning out of control. Eventually the host population got too small to sustain the parasitoids, so the parasitoids went extinct and the host population climbed to its carrying capacity.

Under these parameter conditions, the model is unstable. There are, however, conditions where the model is stable, where neither population goes extinct.

One way we can make it stable is by making \(K\) really small.

P=NULL  #create an empty set for the parasitoids
H=NULL  #create an empty set for the hosts
P[1]=1  #define initial conditions for the parasitoids
H[1]=5  #define initial conditions for the hosts

time=150  #how many time steps to run our simulation
r=0.25  #intrinsic growth rate of the host
a=0.25  #search efficiency of the parasitoid
K=10  #carrying capacity of the host

#here's our model
for(i in 1:time){
  H[i+1]=H[i]*exp(r*(1-H[i]/K)-a*P[i])
  P[i+1]=H[i]*(1-exp(-a*P[i]))
}

#this code below will plot P and H on the y-axis, and time on the x-axis
plot(1:(time+1),H,xlab="generations",ylab="density",ylim=c(0,max(H)))
lines(1:(time+1),H,pch=1)
lines(1:(time+1),P,col="blue")
points(1:(time+1),P,pch=2,col="blue")
legend(.8*time,.9*max(c(H,P)),legend=c("H","P"),col=c("black","blue"),lty=c(1,1),pch=c(1,2))

Notice how now the osciallation get smaller and smaller, eventually moving toward some equilibrium. While there are mathematical solutions to find the equilibrium value, we can estimate it by running the model out many more generations and seeing what number H and P settle on.

Give it a try, and see if there are other ways you can make the model stable.

I hope this was helpful!


Useful references:

Hassell, M.P., 1978. The Dynamics of Arthropod Predator-prey Systems. Princeton University Press.

Murdoch, W.W., Briggs, C.J., Nisbet, R.M., 2003. Consumer-resource Dynamics. Princeton University Press.


This lesson was written by Joe M Kaser using Rmarkdown