Discrete Dynamical Modeling of Ecological Systems

Alexander Jenson

require(pracma)

Part I: Northern Spotted Owl Population

We model the population dyamics of the spotted owls using a Leslie Matrix (named after biologist Patrick Leslie). The population is partitioned into groups based on age classes. The age classes for the spotted owls are:

  • juvenile (less than 1 year old)
  • subadult (1 to 2 years old)
  • adult (over 2 years old)

We assume that the population is evenly split between males and females. This allows us to only track the female population.

The female population at time \(t\) is described by the vector \({\bf x}_t\), given by \[{\bf x}_t=\begin{bmatrix} j_t \\ s_t \\ a_t \end{bmatrix},\]

where \(j_t\) is the number of juveniles, \(s_t\) is the number of subadults and \(a_t\) is the number of adults. The Leslie Matrix Model (also known as the stage-matrix model) tracks (1) the fecundity (reproduction rate) of reproducing adults and (2) the survival rates of the population groups. To find the population at time \(t+1\), we simply multiply by the Leslie Matrix \(L\):

\[{\bf x}_{t+1}=L{\bf x}_t.\]

Here is the Leslie Matrix for owl populations given in the book: \[L=\begin{bmatrix} 0&0&0.33 \\ 0.18 & 0 & 0 \\ 0 & 0.71 & 0.94 \end{bmatrix}.\]

Each time interval represents one year in this case. The first row of the Leslie matrix captures the fecundity of mature spotted owls. This row says that the birth rate for an adult female is 0.33 juveniles per year. In other words, on avergage, an adult female gives birth to a juvenile once every three years.

The second row captures the transition from juvenile to subadult. This row says that each year, 18% of juveniles survive and become subadults.

The final row captures the new adult population. We see that 71% of subadults survive and become adults. Meanwhile, 94% of adults survive the year, and remain in the adult population.

Population forecast

Suppose that we started out with 200 juveniles, 300 subadults and 500 adults. What happens to this population over time?

First, here is a function that computes the trajectory from an initial population distribution. The input variable history can be a single \(3 \times 1\) vector or a \(3 \times T\) matrix with a population history to date, in which case each column represents the data for one year. The variable n denotes how many years more you want to compute in your trajectory.

run.n.steps=function(A,history,n){
  history.dim=dim(history)
  start.length=history.dim[2]
  new.history=matrix(NA,history.dim[1],start.length+n)
  new.history[,1:start.length]=history
  for (i in 1:n){
    new.history[,start.length+i]=A%*%new.history[,start.length+i-1]
  }
  return(new.history)
}

Here is a second function that plots the trajectories output from the run.n.steps function.

plot.trajectory=function(history){
  history.dim=dim(history)
  matplot(t(history),type="l",lwd=3,col=1:history.dim[1],lty=rep(1,history.dim[1]),ylab="population",xlab="time")
  }
  1. We begin with a matrix \(L\) with the values given above for the Leslie matrix.
L = cbind(c(0,0.18,0), c(0,0,0.71), c(0.33,0,0.94))
  1. Here is the initial population distribution in a \(3 \times 1\) matrix:
initial=as.matrix(c(200,300,500))

Here we compute the populations of each age group over the next 5 years.

output = run.n.steps(L,initial,5)
output
##      [,1] [,2]   [,3]     [,4]      [,5]      [,6]
## [1,]  200  165 225.39 220.3014 214.04203 210.70510
## [2,]  300   36  29.70  40.5702  39.65425  38.52756
## [3,]  500  683 667.58 648.6122 638.50031 628.34481
plot.trajectory(output)

Now we do it for 100 years:

output2 = run.n.steps(L,initial,100)
plot.trajectory(output2)

The populations are all decreasing.

Using eigenvalues and eigenvectors, we can explain the asymptotic behavior of the system.

eigen(L)
## $values
## [1]  0.9835927+0.0000000i -0.0217964+0.2059185i -0.0217964-0.2059185i
## 
## $vectors
##               [,1]                  [,2]                  [,3]
## [1,] 0.31754239+0i  0.6820937+0.0000000i  0.6820937+0.0000000i
## [2,] 0.05811107+0i -0.0624124-0.5896338i -0.0624124+0.5896338i
## [3,] 0.94646180+0i -0.0450520+0.4256233i -0.0450520-0.4256233i
lambda1=0.9835927
lambda2=-0.0217964+0.2059185i
lambda3=-0.0217964-0.2059185i
eigenvector = cbind(c(0.31754239, 0.05811107, 0.94646180))

Because lamda1 (the only real eigenvalue we got) is less than one, we can assume that all owl populations will decrease (approach zero) over time.

Save the owls!

Conservation biologists would like to save the spotted owls by increasing the survival rate of junveniles. Currently, that value is 0.18, which means that only 18% of juveniles survive to become subadults.

Therefore we must find the value for juvenile survival rate (with 4 decimal precision) that leads to a stable population.

We plot the population trajectory for \(n=100\), and find the limiting populations of juveniles, subadults and adults.

newL = cbind(c(0,0.2561,0), c(0,0,0.71), c(0.33,0,0.94))
output3 = run.n.steps(newL,initial,100)
plot.trajectory(output3)

The juvenile population approaches 225, the subadult population approaches 88, and the adult population approaches 685.

f=function(x){
  det(cbind(c(0,x,0),c(0,0,0.71),c(0.33,0,0.94))-diag(3))
}
fzero(f,.2)
## $x
## [1] 0.2560819
## 
## $fval
## [1] 0

The function puts the variable ‘x’ where the juvenile survival rate went in matrix L, and too the determinant with a lamda value of 1 (a stable population). This returned a value of 0.2561 for x.

Part II: Blue Whale Population Dynamics

  • The Leslie model considers only the females in the population

  • Time is measured in years (in other settings the time interval may be days or seconds or decades)

  • Individuals are grouped into \(m\) age classes denoted \(C_1,C_2,\ldots,C_m\)

  • An individual’s chance of surviving from one age group to the next depends on its current age group

  • The survival rate \(S_i\) of each age group is known. For \(i \leq m-1\), the value \(S_i\) is the proportion that survive from \(C_i\) to \(C_{i+1}\). Meanwhile, \(S_m\) is the proportion of fully mature adults that survive until the next year

  • The reproduction (fecundity) rate \(F_i\) for each age group is known: \(F_i\) is the average number of offspring born to an individual in \(C_i\)

  • The initial age distribution is known

  • The Leslie matrix for this model is the \(m \times m\) matrix

\[ L= \begin{bmatrix} F_1 & F_2 & F_3 & \ldots & F_{m-1} & F_m \\ S_1 & 0 & 0 &\ldots & 0 &0 \\ 0 & S_2 & 0 &\ldots & 0 &0 \\ \vdots& \ddots& \ddots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & S_{m-2} & 0 & 0 \\ 0 & 0& 0&0 & S_{m-1} & S_m \end{bmatrix}\]

  • It has the vector of fecundities (\(F_1\), . . . , \(F_m\)) across the first row, and it has the vector of survival rates (\(S_1,\ldots , S_{m-1}\)) along the diagonal just below the main diagonal. Finally, the entry in the rightmost corner is \(S_m\)

Whale dynamics

In the 1930’s (before its virtual extinction and a great change in its survival rates), a researhcer studied the blue whale population. Due to the long gestation period, mating habits and migration of the blue whale, a female can produce a calf only once in a two-year period. Thus the age classes for the whale were assumed to be:

  • 0-2 years
  • 2-3 years
  • 4-5 years
  • 6-7 years
  • 8-9 years
  • 10-11 years
  • 12 or more years
  1. We construct a Leslie matrix \(B\) with the following parameters:
  • The vector of fecundities is \([0,0,0.19,0.44,0.5,0.5,0.45]\)
  • The survival rate for all whales 11 years old or younger is 77%, and for all whales 12 or more years, it is 78%
whales = rbind(c(0,0,0.19,0.44,0.5,0.5,0.45),c(0.77,0,0,0,0,0,0),c(0,0.77,0,0,0,0,0),c(0,0,0.77,0,0,0,0),c(0,0,0,0.77,0,0,0),c(0,0,0,0,0.77,0,0),c(0,0,0,0,0,0.77,0.78))
  1. According to this model, what is happening to the whale population? Show some tractories and explain the overall behavior using eigenvalues and eigenvectors. What is the overall growth rate? If after a long time I choose 100 blue whales randomly, how many am I likely to see from each age group? Relate your answers to the eigenvalues and eigenvectors of \(B\).
l1=1.0072102
l2=-0.4670494
l3=0.3793408
vectors=eigen(whales)
v1 = as.matrix(vectors$vectors[,1])
v2 = as.matrix(vectors$vectors[,4])
v3 = as.matrix(vectors$vectors[,7])
output4=run.n.steps(whales,v1,100)
output5=run.n.steps(whales,v2,100)
output6=run.n.steps(whales,v3,100)
plot.trajectory(output4)

plot.trajectory(output5)

plot.trajectory(output6)

We got dramatically different trajectories for each unique real initial population/eigenvector. For the first, we saw that the whale populations both increase and accelerate over time, at a rate of approximately 0.7% per unit time. Our second and third initial populaitons had all real values, but began with negative populations, and ended up converging at a population zero. We know the population after 100 units time to be:

finalpop=as.matrix(output4[,100])
finalpop
##              [,1]
## [1,] 1.1578638+0i
## [2,] 0.8851729+0i
## [3,] 0.6767039+0i
## [4,] 0.5173320+0i
## [5,] 0.3954940+0i
## [6,] 0.3023504+0i
## [7,] 1.0246450+0i

…so the relative proportion of each age group would look like:

finalpop/sum(finalpop)
##               [,1]
## [1,] 0.23346090+0i
## [2,] 0.17847803+0i
## [3,] 0.13644429+0i
## [4,] 0.10431001+0i
## [5,] 0.07974374+0i
## [6,] 0.06096312+0i
## [7,] 0.20659990+0i

Sustainable harvesting of blue whales

If a population is increasing, one might be interested in harvesting a portion of the population for some purpose. We want this harvest rate to be sustainable, meaning that the total population remains at a constant level. For simplicity, we will not distinguish between classes in the harvest (but this could be added with a minor modification).

If a fraction \(h\) of the population is harvested each year, where \(0 < h < 1\), then the population model becomes

\[{\bf x}_{t+1} = (1-h) B {\bf x}_t.\]

To find a sustainable harvest rate \(h\) for blue whales, we plot a trajectory of the system dynamics that includes the sustainable harvesting in order to confirm that the blue whale population reaches a stable total.

We know that the population of whales increases 0.72102% for each unit time because of our eigenvalue. To make up for this, h would also be 0.72102, so as to make \[{\bf x}_{t+1} = {\bf x}_{t}\]. We then create a matrix newB which is our original matrix times 1-h

newB = (1-0.0072102)*whales
output8=output4=run.n.steps(newB,v1,100)
plot.trajectory(output8)