Bifurcation diagram review

We talked about how to create bifurcation diagrams in lab and lecture. Remember that this is a figure where the stable equilibrium population size (\(N^*\)) is plotted as a function of some parameter. In this case, let’s make a bifurcation diagram for the logistic growth model of population growth, with stable equilibrium population size as a function of carrying capacity \(K\).

Let’s make a bifurcation diagram of \(N^*\) as a function of \(K\) for the logistic growth model. Recall that the logistic growth model looks like this: \[ \frac{dN}{dt} = rN[1-\frac{N}{K}] \]

With a bifurcation diagram, you’re asking: for a given value for my parameter of interest, what would my equiliibrium population size be? In this case, if \(K = 100\), then \(N^* = 100\). If \(K = 1000\), then \(N^* = 1000\), and so on.

Here’s an example of how that works. First, we’ll assign some value to \(r\), the per capita growth rate.

r <- 0.8

Then, we’ll create a vector of different values for \(K\). Remember that you’re interested in what happens to \(N^*\) as a function of \(K\) - so, Kset will be the x-axis of your bifurcation diagram.

Kset <- seq(from = 1, to = 100, by = 1)
head(Kset)
## [1] 1 2 3 4 5 6

Remember that we’ve already set our expectations for what \(N^*\) should be for each value of \(K\). So, for the first six values of Kset that you see from head(Kset), what do you think \(N^*\) should be?

We can try this with a simple for() loop:

# create a vector of time steps
test.tmpts <- seq(from = 0.1, to = 10, length = 10000)

# create a holding vector for the population values
test.simu <- rep(NaN, times = length(test.tmpts))

# store the starting conditions in the holding vector
test.simu[1] <- 1

# store the first value of Kset as K.1
K.1 <- Kset[1]

for(i in 2:length(test.tmpts)) {
  # calculate change in time
  deltat <- test.tmpts[i] - test.tmpts[i-1]
  
  # calculate change in population size
  deltaN <- test.simu[i-1] * r * (1 - test.simu[i-1] / K.1) * deltat 
  
  # calculate the total population size and store it in the holding vector
  test.simu[i] <- deltaN + test.simu[i-1]
}

# print the last value in the test.simu vector (the N*)
test.simu[length(test.simu)]
## [1] 1

So, if our \(K = 1\), then our \(N^* = 1\). You could do this for each value of \(K\) in Kset individually: store as K, run the for() loop, pull out the last value in the holding vector, store that as \(N^*\), or… you could write a nested for() loop. You’ll basically tell R:
1. Create a vector of different values for your parameter (x-axis of the bifurcation diagram)
2. Calculate the N* for each parameter - do this by calculating population size, then pull out the last value of N (ASSUMPTION: N will reach equilibrium at the end of the simulation)
3. Store the N* (y-axis of the bifurcation diagram)

We’ve already created a vector of different values for \(K\) called Kset - now, we have to calculate \(N^*\) for each value of K. We’ll do that using the loop below:

# create a holding vector for N*
Nstarset <- rep(NaN, times = length(Kset))

# create a sequence of timepoints at which to calculate population size (in the inner for() loop)
timepoints <- seq(from = 1, to = 25, length = 1000)

### nested for() loop structure:
#### outer loop: choose a value of K from Kset
##### inner loop: calculate population size at every time step
#### outer loop: pull out the last population size (N*) and store it in Nstarset

# this is the start of the outer loop!
# for each element j from the first to last in Kset
for(j in 1:length(Kset)) {
  
  # assign K the value from Kset for this iteration
  K <- Kset[j]
  
  # create a holding vector for the inner loop for population sizes
  N.simu <- rep(NaN, times = length(timepoints)) 
  # set your initial condition
  N.simu[1] <- 10
  
  # inner loop: calculate population size (N)
  # for each element i from the second to the last in timepoints
  for(i in 2:length(timepoints)){ 
    # store previous population size as N
    N <- N.simu[i-1]
    
    # calculate change in time
    deltat <- timepoints[i] - timepoints[i-1] 

    # calculating the change in population size
    deltaN <- N * r * (1 - N / K) * deltat 
    
    # calculating total population size
    N.simu[i] = N + deltaN
  }

  # pull out the last value from N.simu and put it in Nstarset
  Nstarset[j] <- N.simu[length(timepoints)] 
  
}

Great! Let’s check to make sure the loop did what we wanted it to. We’ll look at Nstarset and make sure it’s filled in with numbers.

head(Nstarset)
## [1] 1 2 3 4 5 6
tail(Nstarset)
## [1]  95  96  97  98  99 100

Great! Looks like it worked. Now we can make a plot for our diagram, where Kset is on the x-axis, and Nstarset is on the y-axis:

plot(x = Kset, y = Nstarset, 
     type = 'l', 
     lwd = 2, 
     xlab = "Carrying Capacity, K", 
     ylab = "Stable Equilibria of Popn Size, N*")

Dr. Moeller’s task to you in lecture was to make a bifurcation diagram for \(r\). What do you think that would look like?

# store value of 100 for K
K <- 100

# vector of values of r to iterate over
rset <- seq(from = 0.1, to = 1, by = 0.01)

# create a holding vector for N*
Nstarset.r <- rep(NaN, times = length(rset))

# create a sequence of timepoints at which to calculate population size (in the inner for() loop)
timepoints <- seq(from = 0.1, to = 1000, length = 10000)

# nested for() loop
for(j in 1:length(rset)) {

  r <- rset[j]

  N.simu <- rep(NaN, times = length(timepoints)) 
  N.simu[1] <- 10
  
  for(i in 2:length(timepoints)){ 
    N <- N.simu[i-1]
    deltat <- timepoints[i] - timepoints[i-1] 
    deltaN <- N * r * (1 - N / K) * deltat 
    N.simu[i] = N + deltaN
  }
  
  Nstarset.r[j] <- N.simu[length(timepoints)] 
  
}

plot(x = rset, y = Nstarset.r, 
     type = 'l', 
     lwd = 2, 
     xlab = "Growth rate, r", 
     ylab = "Stable Equilibria of Popn Size, N*",
     ylim = c(0, 101))