Density regulation amplifies environmentally induced population fluctuations

Model specification and analytical derivations

Let Y₀ and Y_t denote respectively the initial population size (population size at time 0) and the population size at time t(t ≥ 1). We assume a discrete-time stochastic Gompertz model for the population dynamics (Reddingius, 1971; Royama, 1981; Sibly et al., 2005; Dennis et al., 2006; Mutshinda & O’Hara, 2011; Mutshinda, O’Hara & Woiwod, 2009; Mutshinda, O’Hara & Woiwod, 2011; Mutshinda et al., 2019) so that the population size at time t ≥ 1 is given by (1)${\displaystyle Y_t=Y_{t-1}exp\left\{r\left(1-\frac{ln\left(Y_{t-1}\right)}{k}\right)+{\varepsilon }_t\right\}}$ where r > 0 is the intrinsic growth rate (i.e., e^r is the multiplicative population growth rate in the absence of density dependence), whereas the realized growth rate Y_t/Y_t−1 decreases with increasing population size Y_t−1, and k > 0 is the natural logarithm of the carrying capacity or equilibrium population size/density denoted by K. The error term ɛ_t, assumed to be Gaussian with mean zero and variance ${\sigma }_t^2$, typically integrates process stochasticity (demographic and environmental stochasticity) and observation error. The Gompertz model as presented here is phenomenological and only describes population level changes without explicitly accounting for individual level processes or age structure, in contrast with mechanistic models that translate individual parameters into population level consequences.

For the purpose of the present study, we rely on simulated data and so, we assume that the population size is observed without error. In addition, we know the data-generating model when in practice the error term typically involves the effect of model misspecification as well. For small populations, demographic stochasticity may be important, and its inverse scaling with the population size introduces time-dependence in the random error term ${\sigma }_t^2$ (e.g., Saether et al., 2000). Since our focus is on the population dynamics in stationary phase, we assume that the population size is large enough for demographic stochasticity to be unimportant. Therefore, we consider the process error to be entirely due to environmental fluctuations in a stochastic environment with constant variance, so that ${\sigma }_t^2={\sigma }^2$ with the error terms ɛ_t,  t ≥ 1 being serially independent. On the natural logarithmic scale with y_t = ln(Y_t), model (1) becomes (2)${\displaystyle y_t=y_{t-1}+r\left(1-\frac{y_{t-1}}{k}\right)+{\varepsilon }_t.}$

Letting β = 1 − rk⁻¹, we can re-write model (2) as (3)${\displaystyle y_t=r+\beta y_{t-1}+{\varepsilon }_t}$ which is a first-order autoregressive [AR(1)] model with autoregressive coefficient β.

From the expression of β in terms of r and k, it follows that $k=r{\left(1-\beta \right)}^{-1}$, which establishes a one-to-one correspondence between the parameters of the Gompertz model in AR(1) form (Eq. 3) and those of the standard Gompertz model (Eq. 1). However, the AR(1) form is more convenient for analyzing the model’s dynamic behavior as we do next.

We start by analyzing the dynamic behavior of the deterministic Gompertz model in AR(1) form (i.e., Eq. 3 without the error term ɛ_t). We discuss the existence of an equilibrium point and the pace at which the population size/density approaches the equilibrium. We then consider the stochastic version of the model (Eq. 3) for which the equilibrium is a non-degenerate distribution (the stationary distribution) rather than being a single point as in the deterministic model. We analyze the temporal fluctuations of the population size around the mean value of its stationary distribution.

Dynamic behavior of the deterministic Gompertz model

The deterministic part of the Gompertz model in AR(1) form is given by the first-order difference equation (4)${\displaystyle y_t=r+\beta y_{t-1}.}$

If $\left|\beta \right|<1$, the system will eventually converge to an equilibrium state y_∞ at which the population size remains constant, meaning that y_∞ = r + βy_∞. Solving this equation for y_∞ yields the following expression for the non-trivial equilibrium (5)${\displaystyle y_{\infty }=r{\left(1-\beta \right)}^{-1}=k.}$

We can use the relation r = y_∞(1 − β) resulting from (5) to describe the population dynamics in terms of y_∞ as $y_1=y_{\infty }+\beta \left(y_0-y_{\infty }\right)$, y₂ = y_∞ + β²(y₀ − y_∞), and by induction we obtain the following closed-form solution to the first-order difference Eq. (4). (6)${\displaystyle y_t=y_{\infty }+{\beta }^t\left(y_0-y_{\infty }\right).}$

Consequently, y_t will converge to y_∞ provided that $\left|\beta \right|<1$, with faster convergence when $\left|\beta \right|$ is closer to zero, with y_t approaching the equilibrium state y_∞ either monotonically from above or from below when 0 < β < 1 or through damped oscillations when −1 < β < 0.

Following from Eq. (6), y_t relates to y_t−1 as $y_t=y_{\infty }+\beta \left(y_{t-1}-y_{\infty }\right),$ which can be re-arranged as $y_t=\left(1-\beta \right)y_{\infty }+\beta y_{t-1}$. The relationship $y_t=\left(1-\beta \right)y_{\infty }+\beta y_{t-1}$ indicates that y_t is a linear combination of y_∞ and y_t−1 with the autoregressive coefficient β providing a direct measure of the dependence of y_t on y_t−1. As a result, the autoregressive parameter of the Gompertz model in AR(1) form is often referred to as the strength of density dependence (e.g., Hampton et al., 2013; Ponciano, Taper & Dennis, 2018; Messmer, 2019; Peeters et al., 2022). We also adopt this terminology.

Deterministic models of population dynamics draw on the assumption that vital rates are constant over time, so that a single set of input values uniquely determines the output value. However, an obvious feature of the real world is that the environment varies continually. Environmental stochasticity refers to the variability in demographic rates caused by random variations in environmental conditions, whereas sampling variations in independent outcomes of demographic events (births, death and dispersal) among individuals in a finite population produces a different kind of fluctuations in population dynamics known as demographic stochasticity. While environmental stochasticity affects all populations irrespective of size, demographic stochasticity is generally only relevant in small populations due to its inverse scaling with the population size. We refer to Lande, Engen & Saether (2003) for details on statistical methods for estimating demographic and environmental stochasticity from empirical data. We next consider the dynamic behavior of the stochastic Gompertz model and investigate the interaction between density-dependent effects and environmental noise in driving temporal fluctuations in population abundance/density.

Dynamic behavior of the stochastic Gompertz model

The stochastic Gompertz model in AR(1) form is given by Eq. (3). If $\left|\beta \right|<1$ the population size will eventually reach a stationary distribution (or equilibrium distribution), which is the stochastic version of an equilibrium in the deterministic model. In stationary phase, the population growth rate is zero on average. Since the additive error ɛ_t are assumed to be normally distributed, the stationary distribution is also normal with long-run mean value μ_∞ and variance v_∞ determined by their time-invariance property. Taking expectations of both sides of Eq. (3) and applying the time-invariance property of the long-run mean yields the equation μ_∞ = r + βμ_∞ whose solution for μ_∞ is (7)${\displaystyle {\mathrm{\mu }}_{\infty }=r{\left(1-\beta \right)}^{-1}=k}$

Alternatively, taking variances of both sides of (3) and applying the time-invariance property of the long-term variance yields the equation v_∞ = β²v_∞ + σ² whose solution for v_∞ is (8)${\displaystyle v_{\infty }={\sigma }^2/\left(1-{\beta }^2\right).}$

Equations (7) and (8) indicate that both the mean of the stationary distribution, which is identical to the deterministic stable equilibrium value y_∞ = k (the log-carrying capacity), and the stationary variance (i.e., the variance of population time series in stationary phase) depend on the strength, β, of density regulation. In addition, the stationary variance v_∞ is proportional to the environmental variance, σ², and relates to the density feedback in such a way that stronger density regulation (i.e., $\left|\beta \right|$ values close to 1) induces higher variability.

It is worth emphasizing that stationarity is a key requirement for analyzing density dependence in population time series, and Eqs. (7) and (8) only make sense for stationary time series. Therefore, the sensible range for the autoregressive parameter β is the interval (−1, 1). However, since we are interested in the stationary variance which depends on β only through β², we restrict attention on β values in the unit interval 0 ≤ β < 1. From our model assumptions that r and k are both positive and the expression β = 1 − rk⁻¹ connecting the autoregressive coefficient of the Gompertz model in AR(1) form to the parameters r and k of the standard Gompertz model, it follows that 0 ≤ β < 1 corresponds to r ≤ k.

When β = 0, which corresponds to the case r = k in the standard Gompertz model formulation, the population trajectory follows a Gaussian white noise process shifted at r. This process is stationary with constant variance σ². When β = 1, which arises as a limiting case when k tends to infinity in the standard Gompertz model formulation, the dynamics are density independent, and the population trajectory is a random walk process with drift r. This process is not stationary since $\mathrm{V\; ar}\left(y_t\right)=t{\sigma }^2$ is unbounded and $\mathrm{E}\left(y_t\right)=y_0+rt$ is only constant when r = 0 (i.e., for the random walk without drift), while a key condition for (weak) stationarity is that the mean and variance of the time series are constant, the other requirement being that for any integer h, Cov(y_t, y_t+h) depends only on h and not on t.

For density-dependent dynamics where β is statistically different from zero and $\left|\beta \right|<1$, the population fluctuates around the stationary mean μ_∞ with variance σ²/(1−β²), where σ² represents the environmental variance. Therefore, the portion of the stationary variance due to density regulation, herein denoted by ${\sigma }_{dd}^2$, is given by the difference between the stationary variance σ²/(1 −β²) under density-regulation and its counterpart σ² in the absence of density feedback on population growth rate. That is, (9)${\displaystyle {\sigma }_{dd}^2={\beta }^2{\sigma }^2/\left(1-{\beta }^2\right).}$ Consequently, the proportion of the stationary variance v_∞ due to density regulation is (10)${\displaystyle {\phi }_{dd}={\sigma }_{dd}^2/v_{\infty }={\beta }^2.}$

Interestingly, the proportion of the stationary variance v_∞ due to density regulation is simply the square of the density-dependence effect β. This provides a biological interpretation of the autoregressive coefficient of the AR(1) model when applied to population dynamics, besides being a measure of the strength of density dependence. Equation (9) provides a tool for separating the stationary variance of population time series into a density-independent component and a density dependent counterpart.

A key motivation for using the Gompertz model is that it can be written as an AR(1) model, which is simple with known statistical properties and established inferential procedures (Mutshinda & O’Hara, 2010). Because the AR(1) model is a linear model with Gaussian errors, maximum likelihood estimates of r, β , and σ², which are identical to least square estimates, can be obtained by performing a linear regression of y_t on y_t−1 (t = 1, 2, 3, …) using standard statistical packages. However, confidence intervals from standard statistical packages are not valid because the $y_t^{\prime }s$ are not independent due to the autoregressive structure in the model. Dennis & Taper (1994) discuss bootstrap and jackknife procedures for constructing confidence intervals of the AR parameters. In a Bayesian framework, priors distributions can be defined to appropriately constrain the model parameters to more sensible range of values, thereby reducing posterior uncertainty while guaranteeing parameter identifiability (Mutshinda, 2009a; Mwanza, 2010).

A common objection to the Gompertz model is that the growth rate depends, if applicable, only logarithmically on the population density allegedly inducing weaker density feedback than the closely related Ricker model (Ricker, 1954) where the growth rate at time t depends on Y_t−1, the actual population size (Dennis & Taper, 1994). The Stochastic Ricker model is given by (11)${\displaystyle Y_t=Y_{t-1}exp\left\{r\left(1-\frac{Y_{t-1}}{K}\right)+{\varepsilon }_t\right\},}$ where r is the intrinsic growth rate and K is the carrying capacity or equilibrium population density on the scale of untransformed population size, and ɛ_t are zero-mean random shocks to the population growth rate assumed to be normally distributed and serially independent. The Ricker model is widely used for population dynamics modeling, particularly in fisheries.

In the next section, we describe a simulation study designed to corroborate empirically the pattern of increasing stationary variance with increasing strength of density regulation under the Gompertz model and to investigate whether this pattern carries over to the Ricker model.

Report on the simulation study

Data simulation

Given numerical values of the parameters r, k and σ, an initial log population size y₀, and a routine for generating standard normal random variables, one can recursively generate trajectories y₀,  y₁, y₂, …y_n from the Gompertz model with initial population density set to the log-carrying capacity. We simulated replicated population time series over 300 time steps and discarded the first 200 observations to ensure that the last 100 data points come from the stationary distribution. We tuned the parameters of the data-simulating model to mimic different levels of density regulation and different levels of environmental noise. Since the parameter β only depends on the positive parameters r and k through the ratio r/k , one can generate data with different levels of density-regulation (different β values) by fixing one of the two parameters, typically the carrying capacity, and varying the other (e.g., Greenwell & Ng, 1984). Fixing the carrying capacity to 1 amounts to expressing population density in units of the carrying capacity. In our simulations, we fix the log-carrying capacity to 1, so that β = (1 − r) for 0 < r < 1. We simulated data under four different levels of environmental noise, with environmental variance set to 0.10, 0.15, 0.20 and 0.25, and three different values of r namely, 0.8, 0.6, and 0.4, corresponding to β values 0.2, 0.4, and 0.6, respectively.

We graphically checked whether the empirical stationary variance v_∞ increased with increasing the strength of density feedback as implied by the analytical result ((8)). We computed the contribution ${\sigma }_{dd}^2$ of density regulation to the empirical stationary variance v_∞ at different levels of density regulation as the difference between the empirical stationary variance of simulated population trajectories and the assumed environmental variance.

In practice, the model parameters are unknown, and one has to rely on estimates obtained from the model fitting to data. We examined the relationship between the strength of density regulation and the proportion of stationary variance due to density regulation using environmental variance estimates over 100 replicated population trajectories at each of the three levels of density regulation under consideration.

A question that springs to mind is whether the pattern of increasing stationary variance with increasing strength of density feedback inferred under the Gompertz model holds for the Ricker model as well. Since under the Ricker model, unlike the Gompertz model, we do not have an expression separating the stationary variance into contributions from environmental stochasticity and density regulation, we rely on simulations to tackle this question. We fit, with a Bayesian approach (Gelman et al., 2013; Mutshinda et al., 2022), the stochastic Gompertz and stochastic Ricker models to data replicates simulated from the stochastic Gompertz model with different levels of density regulation. For the Gompertz model, we independently assigned a Gamma(1, 1) prior on the intrinsic growth rate r, a standard normal prior independently on the autoregressive coefficient β and an InvGamma(0.1, 0.1) prior on the environmental variance σ². For the Ricker model, we independently assigned a Gamma(1, 1) prior on the on the intrinsic growth rate r and Gamma(0.1, 0.1) priors on the carrying capacity K and the environmental variance σ². We used Markov chain Monte Carlo methods (Gilks, Richardson & Spiegelhalter, 1996) to simulate, via the Bayesian freeware OpenBUGS (Thomas et al., 2006; Mutshinda, 2009b), from the joint posterior distributions. We primarily ran the models in OpenBUGS to assess their convergence both informally through visual inspection of traceplots and autocorrelation plots in OpenBUGS, and formally by looking at the Gelman, Rubin statistics via the R package CODA (Plummer et al., 2006), and found that at least 2,000 iterations were required for convergence. Therefore, we used a burn-in period of 4,000 iterations for both models.

We estimated the contribution ${\sigma }_{dd}^2$ of density dependence to the stationary variance v_∞ by the difference between the empirical stationary variance of simulated population time series and the posterior estimate of the environmental variance, and subsequently derived the proportion φ_dd of stationary variance due to density regulation as ${\phi }_{dd}={\sigma }_{dd}^2/v_{\infty }$.