# Models of population growth in continuous/discrete time

`population-growth.Rmd`

Modeling the dynamics of individual species’ populations is a
longstanding challenge in ecology. `ecoevoapps`

includes
functions to simulate various models of population dynamics, including
models in continuous vs. discrete time, models of exponential
vs. logistic growth, and growth of age-structured populations.

Fundamentally, a population can grow due to births or when individuals immigrate from outside sources, and a population can shrink due to death or when individuals leave (emigrate) away from the population. But when we dig a bit deeper, there are actually lots of nuances for us to consider. For example, do births in our focal species happen in discrete time steps (e.g. all births happen during a breeding season), or is it a continuous process (e.g. growth of a bacterial colony in a flask of nutrient-rich broth)? Do all individuals in the population have an equal chance of giving birth and of dying, or do the rates of reproduction and mortality differ by age or life stage?

The package implements various models that differ their biological interpretations and assumptions. Most of these models focus on “closed” populations, which are populations whose size only grow or shrink due to births and deaths – immigration and emigration are not relevant. The package also includes a meta-population model in which individuals can move from high-quality “source” habitats to poorer quality “sink” habitats; this model is described in detail in a different vignette.

Let’s start with models of single population growth in continuous time.

## Continuous time models

### Exponential growth

This model focuses on a closed population where the birth and death-rates do not depend on population density. This setting results in an exponentially growing population, since there is no check on population growth: the birth and death rates remain constant no matter the size of the population. Given that \(b\) is the per-capita birth rate, and \(d\) is the per-capita death rate, the change in population size \(N\) is as follows:

\[\frac{dN}{dT} = bN - dN = (b-d)N\]

Given that \((b-d)\) is a constant that can be expressed simply as \(r\), we can write the exponential growth equation as follows:

\[\frac{dN}{dT} = rN\]

We can simulate the dynamics of the model using the function
`run_exponential_model()`

```
time_exp <- seq(0, 50, 0.1) # length of time for model simulation
init_exp <- c(N1 = 1) # initial population size
params_exp <- c(r = 0.1) # population growth rate
sim_df_exp <- run_exponential_model(time = time_exp, init = init_exp, params = params_exp)
```

This creates a data frame with columns for time and population size:

```
head(sim_df_exp)
#> time N1
#> [1,] 0.0 1.000000
#> [2,] 0.1 1.010050
#> [3,] 0.2 1.020202
#> [4,] 0.3 1.030455
#> [5,] 0.4 1.040811
#> [6,] 0.5 1.051271
```

We can plot the size of this population over time using the function
`plot_continuous_population_growth()`

:

`plot_continuous_population_growth(sim_df_exp)`

### Logistic growth

The exponential growth model above is very straightforward, and predicts that in the long-term, any population experiencing such growth will reach an infinitely large size. Given that populations never grow without limits in nature, we can improve on this model.

Let us re-examine the assumptions of the exponential growth model– namely, that there is a per-capita birth rate \(b\) and a per-capita death rate \(d\) that stays constant at all times. We can simply relax this assumption by examining what happens when the net growth rate (\(rD\)) declines linearly as a population reaches an environmentally-determined carrying capacity (\(K\)). The dynamics of such a population can be described by the logistic growth equation:

\[\frac{dN}{dt} = rN \left(1-\frac{N}{K}\right)\]

Populations experience logistic growth can grow (i.e., have a positive population growth rate) until the population size is equal to the carrying capacity (\(N = K\)). When a population exceeds its carrying capacity, the population will have a negative growth rate, and shrink back to the carrying capacity. For more on this model, refer to this link.

We can simulate the dynamics of a population with logistic growth
using the function `run_logistic_model()`

as follows:

```
time_log <- seq(0, 100, 0.1) # length of time for model simulation
init_log <- c(N1 = 1) # initial population size
params_log <- c(r = 0.2, K = 100) # population growth rate
sim_df_log <- run_logistic_model(time = time_log, init = init_log, params = params_log)
```

As above, this creates a data frame with columns for time and population size:

```
head(sim_df_log)
#> time N1
#> [1,] 0.0 1.000000
#> [2,] 0.1 1.019996
#> [3,] 0.2 1.040387
#> [4,] 0.3 1.061182
#> [5,] 0.4 1.082387
#> [6,] 0.5 1.104012
```

We can plot the trajectory of this population over time using the same code above:

`plot_continuous_population_growth(sim_df_log)`

### Lagged logistic growth

In some circumstances, the growth rate of a population might depend on the population size at some time in the past. For example, the reproductive output of individuals may be limited by how much they had to compete for resources as juveniles. We can model such a scenario by making the population growth rate at time \(t\) a function of some time \(t-\tau\) in the past as follows:

\[ \frac{dN}{dt} = rN_{t} \left(1-\frac{N_{t-\tau}}{K}\right)\]

We can model the dynamics of a population with lagged logistic growth
using the same `run_logistic_model()`

function as above, but
we now simply need to define the additional parameter `tau`

in the parameters vector:

```
time_laglog <- seq(0, 100, 0.1) # length of time for model simulation
init_laglog <- c(N1 = 1) # initial population size
params_laglog <- c(r = 0.5, K = 100, tau = 2) # population growth rate
sim_df_laglog <- run_logistic_model(time = time_laglog, init = init_laglog, params = params_laglog)
```

This again produces a data frame with columns time and N1, though we are not showing it here to avoid repetition. At sufficiently large values of \(r\) and/or \(\tau\), this model results in damped oscillations in population sizes, before the population ultimately converges to its carrying capacity:

`plot_continuous_population_growth(sim_df_laglog)`

## Discrete time models of population growth

The models above apply well to populations where births and deaths happen continuously, but they are not very relevant to species where births and deaths happen at discrete times (e.g. during a single breeding season). To model such species, we can turn to models of population growth in discrete time.

### Exponential growth in discrete time

The simplest model of population growth in discrete time assumes that the population size at time \(t+1\) (\(N_{t+1}\)) is a product of the population size at time \(t\) (\(N_t\)) and the population growth rate, symbolized \(\lambda\):

\(N_{t+1} = \lambda N_t\)

When \(\lambda > 1\), the
population grows every year, resulting in exponential growth, and when
\(\lambda < 1\), the population
shrinks to extinction. We can simulate the dynamics of such a population
with the function `run_discrete_exponential_model()`

:

```
# We simply provide initial population size, lambda, and the time
# for simulation to the model as arguments:
sim_df_dexp <- run_discrete_exponential_model(N0 = 1, lambda = 1.2, time = 20)
```

This produces a data frame with columns time and N1:

```
head(sim_df_dexp)
#> time Nt
#> 1 0 1.00000
#> 2 1 1.20000
#> 3 2 1.44000
#> 4 3 1.72800
#> 5 4 2.07360
#> 6 5 2.48832
```

We can plot the population growth over time using the function
`plot_discrete_population_growth()`

:

`plot_discrete_population_growth(sim_df_dexp)`

### Logistic growth in discrete time

There are several ways we can build in limits to population growth in discrete time models. One of the most ‘obvious’ ways might be to convert the continuous time logistic growth model we met above into a model of discrete time population growth:

\[N_{t+1} = r_dN_t\left(1-\frac{N_t}{K}\right)\]

In this case, \(r_d\) takes the place of \(\lambda\), because it is the discrete version of the \(r\) parameter in the logistic growth model. Even though this seems like a fairly inocuous switch from continuous to discrete time, it turns out that this model is more interesting than it might look at first sight. In fact, this formulation was notably studied by Robert May in the classic 1976 paper “Simple mathematical models with very complicated dynamics”, in which he showed the potential for this model to generate chaotic dynamics. Another interesting aspect of this model is that unlike in the contiuous time model, \(K\) here doesn’t represent “carrying capacity” in the traditional sense, since the population doesn’t actually equilibrate to that population size – it equilibrates to a population size lower than \(K\). For more on this model, refer to this link.

We can simulate the model using the function
`run_discrete_logistic_model()`

, as follows:

```
params_disclog <- c(rd = 1.5, K = 1000)
sim_df_disclog <- run_discrete_logistic_model(N0 = 1, params = params_disclog, time = 50)
```

As above we can use the
`plot_discrete_population_growth()`

function to plot the
dynamics of this model:

`plot_discrete_population_growth(sim_df_disclog)`

Robert May’s 1974 paper (show above) showed that at sufficiently
large values of `$r_d$`

, this simple model was capable of
producing not only cyclic population dynamics, but also chaotic
dynamics:

```
# By default, the run_discrete_logistic_model uses N0 = 1, so we can omit it here:
params_disclog_cyc = c(rd = 3, K = 1000)
sim_df_disclog_cyc <- run_discrete_logistic_model(params = params_disclog_cyc, time = 75)
params_disclog_chaos <- c(rd = 3.75, K = 1000)
sim_df_disclog_chaos <- run_discrete_logistic_model(params = params_disclog_chaos, time = 75)
plot_discrete_population_growth(sim_df_disclog_cyc) +
plot_discrete_population_growth(sim_df_disclog_chaos)
```

One helpful way to visualize the dynamics of a discrete logistic model is with a “cobweb plot”, which shows how the population at time \(t\) relates to the population at \(t+1\). For more information on cobweb plots, how to read them, and what they tell us, here are some useful references: link 1; link 2; link 3.

We can make a cobweb plot using the function
`plot_discrete_population_cobweb()`

. For this function, we
need to provide the simulated data frame, as well as the parameter
vector and the model type:

```
plot_discrete_population_cobweb(sim_df_disclog_cyc,
params_vec = params_disclog_cyc,
model_type = "discrete_logistic")
#> Warning: Removed 1 rows containing missing values (geom_segment).
#> Removed 1 rows containing missing values (geom_segment).
```

Compare the cobweb plot of the population with cyclic dynamics with that of the population with chaotic behavior:

```
plot_discrete_population_cobweb(sim_df_disclog_chaos,
params_vec = params_disclog_chaos,
model_type = "discrete_logistic")
#> Warning: Removed 1 rows containing missing values (geom_segment).
#> Removed 1 rows containing missing values (geom_segment).
#> Warning: Removed 15 row(s) containing missing values (geom_path).
```

### Other models of logistic growth in discrete time

As suggested above, there are several other commonly used models of
population growth in discrete time, each with distinct behaviors and
assumptions. Two models that are supported in `ecoevoapps`

are the Beverton Holt model and the Ricker model. Both of these were
developed in the context of fisheries management but have since been
applied to various other organisms.

#### Beverton Holt model

The Beverton-Holt model models population dynamics as follows:

\[N_{t+1} = \frac{RN_t}{1+\left(\frac{R-1}{K}\right)N_t}\]

#### Ricker model

The Ricker model models population dynamics as follows: \[N_{t+1} = N_t e^{(r (1 - N_t/K))}\]

#### Simulating in `ecoevoapps`

The Beverton-Holt model and Ricker model can be simulated in
`ecoevoapps`

using the functions
`run_beverton_holt_model()`

and
`run_ricker_model()`

respectively. The usage is same as in
the discrete logistic model above, so we won’t repeat the code here for
brevity. To make cobweb plots, note that the parameter
`model_type`

in `plot_discrete_population_cobweb`

needs to be set to `beverton_holt`

or `ricker`

as
appropriate.

## Structured population growth

The models above assume that within populations, all individuals have
equal rates of reproduction and death. But in fact, some populations are
best modeled as a collection of individuals in distinct age classes (Tenhumberg,
2010). The growth of such structured populations depends on the
transition rates from one age class to another. For example, the
survival rates of juveniles into young adults may be different from the
survival rate of young adults into older adults. We can model the
dynamics of such populations by knowing the rates of *survival*
from one age class into the next, and by knowing how individuals in each
age contribute to the population by giving birth to new individuals.

### The Leslie matrix

The Leslie
Matrix is a convenient way of organizing age transition rates
(survival and fecundity). The Leslie matrix is an \(NxN\) matrix, where \(N\) indicates the total number of ages in
the population cycle. The elements of the matrix tell us about how the
individuals in each age *column* contribute to each age
*row*. For example, the second column in the first row tells us
how individuals who are of age 2 contribute to the age 1 category.

In addition to being a convenient way to account for survival and fecundity information, the matrix also has three key properties:

To simulate the population at the next time step (time \(t+1\)), one simply needs to multiply the Leslie matrix by the matrix with the initial population sizes (at time \(t\)).

The largest eigenvalue of the Leslie matrix is \(\lambda\), the asymptotic growth rate for the whole population. When \(\lambda < 1\) the declines to extinction, \(\lambda > 1\) the grows exponentially, and when \(\lambda = 1\) the population size remains constant through time).

- The eigenvector that corresponds to the largest eigenvalue gives the stable population structure, i.e., the proportions of each age class after the population reached equilibrium.

### Simulating structured population dynamics in R

We can model the dynamics of an age-structured population using the
function `run_structured_population_simulation()`

as
follows:

```
# define the Leslie matrix:
leslie_mat <- matrix(c(0, 8, 1, 1,
0.4, 0, 0, 0,
0, 0.8, 0, 0,
0, 0, 0.1, 0),
ncol = 4, byrow = T)
# define a vector of initial population sizes in each class:
init_vec <- c(10, 0, 0, 0)
structured_popgrowth <- run_structured_population_simulation(leslie_mat = leslie_mat, init = init_vec, time = 50)
```

Before we visualize the model output, we can make a population
structure diagram for the defined Leslie matrix using the function
`plot_leslie_diagram()`

:

`plot_leslie_diagram(leslie_mat)`

```
#> $arr
#> row col Angle Value rad ArrowX ArrowY TextX TextY
#> 1 2 1 0 0.4 0.125 0.2494110 0.2750014 0.250 0.255
#> 2 1 2 0 8 0.125 0.2501963 0.5249998 0.250 0.545
#> 3 3 2 0 0.8 0.125 0.4994110 0.2750014 0.500 0.255
#> 4 1 3 0 1 0.250 0.3753927 0.6499997 0.375 0.670
#> 5 4 3 0 0.1 0.125 0.7494110 0.2750014 0.750 0.255
#> 6 1 4 0 1 0.375 0.5005890 0.7749995 0.500 0.795
#>
#> $comp
#> x y
#> [1,] 0.125 0.4
#> [2,] 0.375 0.4
#> [3,] 0.625 0.4
#> [4,] 0.875 0.4
#>
#> $radii
#> x y
#> [1,] 0.1 0.2270617
#> [2,] 0.1 0.2270617
#> [3,] 0.1 0.2270617
#> [4,] 0.1 0.2270617
#>
#> $rect
#> xleft ybot xright ytop
#> [1,] 0.025 0.1729383 0.225 0.6270617
#> [2,] 0.275 0.1729383 0.475 0.6270617
#> [3,] 0.525 0.1729383 0.725 0.6270617
#> [4,] 0.775 0.1729383 0.975 0.6270617
```

This function generates a matrix with as many rows as the number of age classes defined in the leslie matrix, and as many columns as the number of time steps (plus the initial time step):

```
dim(structured_popgrowth)
#> [1] 4 51
```

We can plot the number of individuals in each age class using the
function `plot_structured_population_agedist()`

:

```
plot_structured_population_size(structured_popgrowth)
#> Warning: Transformation introduced infinite values in continuous y-axis
```

We can also plot the proportion of individuals in each age class
using the function
`plot_structured_population_agedist()`

:

`plot_structured_population_agedist(structured_popgrowth)`