Exercise: Structured populations - Solutions

Population projection

Note that this exercise consists of two tabs!

Using the example from the lecture, you will in this exercise apply matrix multiplication to project population dynamics for a couple of time steps into the future.

The following code shows how to define a population projection matrix in R. We here use the three-by-three matrix from the lecture.

transition.matrix <- matrix(c(0.1,0,100,0.02,0.2,0,0,0.2,0.6), nrow=3, ncol=3, byrow=T)
transition.matrix # let's have a look at the matrix
##      [,1] [,2]  [,3]
## [1,] 0.10  0.0 100.0
## [2,] 0.02  0.2   0.0
## [3,] 0.00  0.2   0.6

Let us now define the initial population vector at time \(t=0\).

N.0 <- c(1000, 30, 15)
N.0 # let's have a look at the vector
## [1] 1000   30   15

In order to obtain the vector describing the population composition one time step into the future, we multiply the current population vector N.0 from the right to the population projection matrix transition.matrix. The multiplication between a matrix and a vector is executed in R using the symbol %*%.

##      [,1]
## [1,] 1600
## [2,]   26
## [3,]   15
##        [,1]
## [1,] 1660.0
## [2,]   37.2
## [3,]   14.2

Use the this procedure to produce the population vectors for the first ten time steps, N.1, N.2,..., N.10.

N.3 <- transition.matrix %*% N.2
N.4 <- transition.matrix %*% N.3
N.5 <- transition.matrix %*% N.4
N.6 <- transition.matrix %*% N.5
N.7 <- transition.matrix %*% N.6
N.8 <- transition.matrix %*% N.7
N.9 <- transition.matrix %*% N.8
N.10 <- transition.matrix %*% N.9

You can see the correct results in the following table.

##    time    Seeds      Juv       Ad
## 1     0 1000.000 30.00000 15.00000
## 2     1 1600.000 26.00000 15.00000
## 3     2 1660.000 37.20000 14.20000
## 4     3 1586.000 40.64000 15.96000
## 5     4 1754.600 39.84800 17.70400
## 6     5 1945.860 43.06160 18.59200
## 7     6 2053.786 47.52952 19.76752
## 8     7 2182.131 50.58162 21.36642
## 9     8 2354.855 53.75894 22.93617
## 10    9 2529.103 57.84888 24.51349
## 11   10 2704.259 62.15183 26.27787

Next, let us compute the total number of individuals at time t=0, t=1 until t=10 using the command sum().

Ntot.0 <- sum(N.0) 
Ntot.1 <- sum(N.1)
Ntot.2 <- sum(N.2)
Ntot.3 <- sum(N.3)
Ntot.4 <- sum(N.4)
Ntot.5 <- sum(N.5)
Ntot.6 <- sum(N.6)
Ntot.7 <- sum(N.7)
Ntot.8 <- sum(N.8)
Ntot.9 <- sum(N.9)
Ntot.10 <- sum(N.10)

You can see the correct values in the last column of the following table.

##    time    Seeds      Juv       Ad     Ntot
## 1     0 1000.000 30.00000 15.00000 1045.000
## 2     1 1600.000 26.00000 15.00000 1641.000
## 3     2 1660.000 37.20000 14.20000 1711.400
## 4     3 1586.000 40.64000 15.96000 1642.600
## 5     4 1754.600 39.84800 17.70400 1812.152
## 6     5 1945.860 43.06160 18.59200 2007.514
## 7     6 2053.786 47.52952 19.76752 2121.083
## 8     7 2182.131 50.58162 21.36642 2254.079
## 9     8 2354.855 53.75894 22.93617 2431.550
## 10    9 2529.103 57.84888 24.51349 2611.465
## 11   10 2704.259 62.15183 26.27787 2792.689

Now we have everything to create a time series plot. With this we mean a plot with time on the x-axis and Ntot on the y-axis.

plot(x=c(0:10), y=c(Ntot.0, Ntot.1, Ntot.2, Ntot.3, Ntot.4, Ntot.5, Ntot.6, Ntot.7, Ntot.8, Ntot.9, Ntot.10), type = "b", xlab="Time (years)", ylab="Ntot.t", las=1)

Next, let’s have a look at the population composition. How is the total number of individuals distributed over the different life-history states?

barplot(cbind(N.0, N.1, N.2, N.3, N.4, N.5, N.6, N.7, N.8, N.9, N.10), las=1, col=c("grey", "green", "red"), names.arg=c("N.0", "N.1", "N.2", "N.3", "N.4", "N.5", "N.6", "N.7", "N.8", "N.9", "N.10"), ylab="Number of individuals", xlab="Time (years)", ylim = c(0, 3000), cex.names=0.8)

legend("topleft", bty="n", legend=c("Seeds", "Seedlings", "Adults"), fill=c("grey", "green", "red"), horiz=T)

It is also common (and informative) to plot the population vectors as relative numbers to show what relative proportion of the population belongs to each life history state. For this we have to divide the number of individuals in each state by the total population size at that time step.

barplot(cbind(N.0/sum(N.0), N.1/sum(N.1), N.2/sum(N.2), N.3/sum(N.3), N.4/sum(N.4), N.5/sum(N.5), N.6/sum(N.6), N.7/sum(N.7), N.8/sum(N.8), N.9/sum(N.9), N.10/sum(N.10))*100, col=c("grey", "green", "red"), las=1, names=c("N.0", "N.1", "N.2", "N.3", "N.4", "N.5", "N.6", "N.7", "N.8", "N.9", "N.10"), ylab="Number of individuals", xlab="Time (years)", cex.names=0.8)
legend("bottom", bty="n", legend=c("Seeds", "Seedlings", "Adults"), fill=c("grey", "green", "red"), horiz=T)

Finally, let us calculate the growth factor \(\lambda\) at each time step. For this we have to divide the total population size of two consecutive time steps by each other: \[ \lambda(t) = \frac{N_{\rm tot}(t+1)}{N_{\rm tot}(t)}. \]

You can see the correct values in the last column of the following table.

##    time    Seeds      Juv       Ad     Ntot   lambda
## 1     0 1000.000 30.00000 15.00000 1045.000 1.570335
## 2     1 1600.000 26.00000 15.00000 1641.000 1.042901
## 3     2 1660.000 37.20000 14.20000 1711.400 0.959799
## 4     3 1586.000 40.64000 15.96000 1642.600 1.103222
## 5     4 1754.600 39.84800 17.70400 1812.152 1.107806
## 6     5 1945.860 43.06160 18.59200 2007.514 1.056572
## 7     6 2053.786 47.52952 19.76752 2121.083 1.062702
## 8     7 2182.131 50.58162 21.36642 2254.079 1.078733
## 9     8 2354.855 53.75894 22.93617 2431.550 1.073992
## 10    9 2529.103 57.84888 24.51349 2611.465 1.069395
## 11   10 2704.259 62.15183 26.27787 2792.689       NA

Let’s make a plot that shows how \(\lambda\) changes over time.

plot(x=c(0:9), y=c(lambda.0, lambda.1, lambda.2, lambda.3, lambda.4, lambda.5, lambda.6, lambda.7, lambda.8, lambda.9), type = "b", xlab="Time (years)", ylab="growth factor lambda", las=1)

Let us now repeat the calculations for the two other starting vectors. Note that I here use two other starting vectors than what is stated in the exercise. However, the conslusions are independent of this. First, \(N(0)=(1000, 0, 0)\),

N.0 <- matrix(c(1000, 0, 0), nrow=3, ncol=1, byrow=T)

and then \(N(0)=(0, 30, 0)\).

N.0 <- matrix(c(0, 30, 0), nrow=3, ncol=1, byrow=T)

Answer the following questions.

  1. Project the population ten years into the future and describe in words how the population changes in both absolute. Does the population fluctuate, increase or decrease?

Population size fluctuates until the fourth time step and from then on shows a monotonic increase. If we would project the time series for more than ten time steps, we could see that the increase in total population size is exponential. The distribution of individuals over the three stages also fluctuates slightly during the first time steps but is almost constant from the fourth time step onward. This is in fact an important property of matrix models. Once the population starts to grow (or decline) exponentially the population also reaches a stable stage distribution: the proportion of individuals in the different stages stays constant, here given by (0.968, 0.022, 0.01).

  1. Compute the growth factor \(\lambda\) for ten years and describe how it changes over time. What do you observe.

The observation under (1) is reflected in the behavior of \(\lambda\): this annual growth factor fluctuates slightly for the first four time steps but then stabilizes around 1.07. In other words, from the fourth time step onward, the population grows from one year to the next by a fixed factor 1.07.

  1. Execute the computations under (1) and (2) for three different initial population vectors \(N_0\), namely: \(N_0=(1000, 30, 15)\), \(N_0=(1000, 0, 0)\) and \(N_0=(0, 30, 0)\). Can you make a figure that shows the time series for the three different starting values together? Similarly, can the time series for the three different starting values together in the same figure? Can you do a similar figure but with the \(\lambda\)-values as computed from the different starting vectors? Summarize your observations.

For the alternative starting vectors, we can see that the initial fluctuations in the \(\lambda\)-values are quite extreme but that from time step four onward, the fluctuations disappear and \(\lambda\)-values become very similar. In all cases, \(\lambda\approx 1.07\) after some time. The starting vectors do, however, cause differences in the values for \(N_{\rm tot}\). However, eventually in all three cases, populations will grow exponentially.

Conservation measures

Consider a population of a rare meadow plant. In this exercise, we explore how we can use matrix models to investigate the efficiency of different conservation strategies. Assume that after an intensive period of data collection we determined the following population transition matrix.

meadow.matrix <- matrix(c(0.001,0,200,0.01,0.2,0,0,0.2,0), nrow=3, ncol=3, byrow=T)
meadow.matrix
##       [,1] [,2] [,3]
## [1,] 0.001  0.0  200
## [2,] 0.010  0.2    0
## [3,] 0.000  0.2    0

and the initial population vector N.0 given by

N.0 <- matrix(c(3000, 180, 30), nrow=3, ncol=1, byrow=T)
N.0
##      [,1]
## [1,] 3000
## [2,]  180
## [3,]   30

At the end of this exercise you should answer the following questions.

  1. Draw the life cycle graph corresponding to this matrix. What kind of life history does this plant have?
Life cycle graph

Life cycle graph

This plant is an semelparous, biennal plant with seed bank.

  1. Calculate the projected population growth rate \(\lambda\) and the total population size N.tot for 10 years and plot them as a function of time. After how many years would this population be expected to have less than 10 adults where it would be considered severely endangered?

See graphs below. Note that for this matrix and this initial population vector \(\lambda\) is still somewhat fluctuating after ten time steps. However, after 10 time steps it becomes already clear that the stable value of \(\lambda\) will be below 1 and that therefore without conservation measurements this population will go extinct. According to the below graphs the adult population size drops below 10 after five time steps.

  1. To protect this population, two different measures are considered. For each of these methods, project the population dynamics and plot the number of adults against time. Which conservation measure is preferable?
  1. Adding 100 adult plants from a different population once at year \(t=0\).

This conservation measurement only has a short term benefit but ultimately this population decreases as well. In fact, we could have anticipated this result since we found out in the first part of this exercise that the long term dynamics of a population are independent of the initial population vector and only depend on the entries of the transition matrix.

  1. Increasing recruitment of juveniles from seeds by creating small open sites in the meadow each year. This is expected to result in a transition probability from seeds to juveniles equal to \(0.025\).

This conservation measurement results in a long term increase of the population. Thus, this is clearly the preferred conservation measurement. This must mean that the second conservation measurement caused that \(\lambda\) becomes a number larger than 1.

In order to make investigating these questions easier, we need some additional R-magic. In particular, we want to use R-code that let’s us compute a whole time series with a single command. For this, we make use of what is called a for-loop.

Before we “run” the for-loop to create a time series, we create a so-called data frame in which we can store the data from the time series of length 10. We name this data frame pop.dev, abbreviating population development, and initially fill all entries with NA(for “not available”). You can think of a data frame as an EXCEL sheet but in R. Thus, a data frame is a matrix where we can give names to rows and columns to describe what is stored in the different cells of the data frame.

pop.dev <- data.frame(time=c(0:10), Seeds=rep(NA,11), Juv=rep(NA,11), Adults=rep(NA,11), N=rep(NA,11), lambda=rep(NA,11))

Let us fill the first column of the data frame with the population vector at time zero.

pop.dev[1, "Seeds"] <- N.0[1]
pop.dev[1, "Juv"] <- N.0[2]
pop.dev[1, "Adults"] <- N.0[3]
pop.dev[1, "N"] <- sum(N.0)

Next, let us “wrap” a loop of length 10 around the matrix multiplication.

N.t <-N.0 # this creates a new vector, called N.t, and we give it the values that are stored in N.0
for (n in 2:11) # the loop "runs" from n=2 to n=11 (representing the years 1 to 10 in the data frame)
{
  N.t.plus.1 <- meadow.matrix %*% N.t # computes the population vector for the next time step
  pop.dev[n, "N"] <- sum(N.t.plus.1) # the following lines enter the outcome into the data frame
  pop.dev[n, "Seeds"] <- N.t.plus.1[1]
  pop.dev[n, "Juv"] <- N.t.plus.1[2]
  pop.dev[n, "Adults"] <- N.t.plus.1[3]
  pop.dev[n, "lambda"] <- sum(N.t.plus.1)/sum(N.t)
  N.t <- N.t.plus.1 # updates the current values stored in N.t
}

The “filled” data frame should look as follows.

pop.dev

With the following command, we can plot the numbers for a specific column in the data frame, for example, the total population size N.

plot(pop.dev$time, pop.dev$N, type="l", ylab="Total population size", xlab="Time (years)", las=1)

In order to plot the entries of another column, we have to add the name of that column to the second pop.dev$ in the last line of code. For example, if we want to plot the number of adults we have to replace pop.dev$N with pop.dev$Adults (and change the name of the label on the y-axis accordingly).

plot(pop.dev$time, pop.dev$Adults, type="l", ylab="Number of Adults", xlab="Time (years)", las=1)

And here is the code if you wish to plot the growth factor \(\lambda\) as a function of time.

plot(pop.dev$time, pop.dev$lambda, type="l", ylab="Growth factor lambda", xlab="Time (years)", las=1)

Now you have to decide how you have to change the above code so that it takes into account the above two conservation strategies.

The following code produces the data for the FIRST conservation measurement. Here, I only changed the vector N.0, by adding 100 adult individuals, from \(N_0=(3000, 180, 30)\) to \(N_0=(3000, 180, 130)\).

N.0 <- matrix(c(3000,180,130), nrow=3, ncol=1, byrow=T)

pop.dev.strategy.1 <- data.frame(time=c(0:10), Seeds=rep(NA,11), Juv=rep(NA,11), Adults=rep(NA,11), N=rep(NA,11), lambda=rep(NA,11))

pop.dev.strategy.1[1, "Seeds"] <- N.0[1]
pop.dev.strategy.1[1, "Juv"] <- N.0[2]
pop.dev.strategy.1[1, "Adults"] <- N.0[3]
pop.dev.strategy.1[1, "N"] <- sum(N.0)

N.t <-N.0 # the initial value of N.n is given by the value N.0
for (n in 2:11) # the year runs second to 11th row (i.e, year 1 to 10)
{
  N.t.plus.1 <-  meadow.matrix %*% N.t
  pop.dev.strategy.1[n, "N"] <- sum(N.t.plus.1)
  pop.dev.strategy.1[n, "Seeds"] <- N.t.plus.1[1]
  pop.dev.strategy.1[n, "Juv"] <- N.t.plus.1[2]
  pop.dev.strategy.1[n, "Adults"] <- N.t.plus.1[3]
  pop.dev.strategy.1[n, "lambda"] <- sum(N.t.plus.1)/sum(N.t)
  N.t <- N.t.plus.1
}

plot(pop.dev.strategy.1$time, pop.dev.strategy.1$Adults, type="l", ylab="Number of Adults", xlab="Time (years)", las=1, col="green")

The following code produces the data for the SECOND conservation measurement. For this we have to change the first entry in the second column of population projection matrix from 0.01 to 0.025. This entry describes the probability of seeds to develop into small plants.

meadow.matrix <- matrix(c(0.001,0,200,0.025,0.2,0,0,0.2,0), nrow=3, ncol=3, byrow=T)
N.0 <- matrix(c(3000,180,30), nrow=3, ncol=1, byrow=T)

pop.dev.strategy.2 <- data.frame(time=c(0:10), Seeds=rep(NA,11), Juv=rep(NA,11), Adults=rep(NA,11), N=rep(NA,11), lambda=rep(NA,11))

pop.dev.strategy.2[1, "Seeds"] <- N.0[1]
pop.dev.strategy.2[1, "Juv"] <- N.0[2]
pop.dev.strategy.2[1, "Adults"] <- N.0[3]
pop.dev.strategy.2[1, "N"] <- sum(N.0)

N.t <-N.0 # the initial value of N.n is given by the value N.0
for (n in 2:11) # the year runs second to 11th row (i.e, year 1 to 10)
{
  N.t.plus.1 <- meadow.matrix %*% N.t
  pop.dev.strategy.2[n, "N"] <- sum(N.t.plus.1)
  pop.dev.strategy.2[n, "Seeds"] <- N.t.plus.1[1]
  pop.dev.strategy.2[n, "Juv"] <- N.t.plus.1[2]
  pop.dev.strategy.2[n, "Adults"] <- N.t.plus.1[3]
  pop.dev.strategy.2[n, "lambda"] <- sum(N.t.plus.1)/sum(N.t)
  N.t <- N.t.plus.1
}

plot(pop.dev.strategy.2$time, pop.dev.strategy.2$Adults, type="l", ylab="Number of Adults", xlab="Time (years)", las=1, col="blue")

Once you have the time series for the three different scenarios (population without conservation, population with measurement 1 and population with measurement 2), you could use code like the following to visualize the results in a single figure.

plot(pop.dev$time, pop.dev$Adults, type="l", ylab="Number of Adults", xlab="Time (years)", las=1, col="red", ylim = c(0, 60))
lines(pop.dev$time, pop.dev.strategy.1$Adults, col="green")
lines(pop.dev$time, pop.dev.strategy.2$Adults, col="blue")
legend("top", legend=c("no measure", "Meth. 1", "Meth. 2"), 
       col=c("red", "green", "blue"), lty=c(1,1,1), bty="n", horiz=TRUE, cex=0.9)