The diversity of population responses to environmental change

Abstract The current extinction and climate change crises pressure us to predict population dynamics with ever‐greater accuracy. Although predictions rest on the well‐advanced theory of age‐structured populations, two key issues remain poorly explored. Specifically, how the age‐dependency in demographic rates and the year‐to‐year interactions between survival and fecundity affect stochastic population growth rates. We use inference, simulations and mathematical derivations to explore how environmental perturbations determine population growth rates for populations with different age‐specific demographic rates and when ages are reduced to stages. We find that stage‐ vs. age‐based models can produce markedly divergent stochastic population growth rates. The differences are most pronounced when there are survival‐fecundity‐trade‐offs, which reduce the variance in the population growth rate. Finally, the expected value and variance of the stochastic growth rates of populations with different age‐specific demographic rates can diverge to the extent that, while some populations may thrive, others will inevitably go extinct.


Growth rates for populations with same deterministic λ and average adult rates but different number of age classes
A. Simulation of mortality and fecundity trajectories. We constructed five mortality rates and five fecundity rates that covered 25 combinations of age-specific survival and fecundity trajectories depicted in the left panel of Fig. 3 in the main text.
We calculated the range of survival trajectories by varying the mortality or hazard rate function, defined as where x is age and X is a random variable for ages at death, while θ is a vector of parameters. From Eq. S1a, several demographic functions are derived, namely S(x|θ) = Pr (X > x) = exp − x 0 µ(y|θ)dy , [S1b] F (x|θ) = Pr (X < x) = 1 − S(x|θ), [S1c] where Eq. S1b is the survival probability, Eq. S1c is the the probability that death occurs before age x (or cumulative density function, CDF), and Eq. S1d is the probability density function (PDF) of ages at death. We used five baseline age-specific mortality functions where the first (M1) produced a declining mortality with age (47) given by µ(x|θ) = exp(a0 − a1x) + c, where a0 ∈ R and a1, c > 0. The second (M2) produced a bathtub or Siler mortality profile (39) as where a0, b0 ∈ R and a1, b1, c > 0. The third (M3) produced typical senescent profile with exponentially increasing mortality given by the Gompertz-Makeham model (18,24) given by where b0 ∈ R and b1, c > 0. The forth (M4) and fifth (M5) mortality models were constructed using a logistic mortality model (31) of the form µ(x|θ) = c + exp(b0 − b1x) , where b0 ∈ R and b1, b2, c > 0. The main difference between M4 and M5 was given by the intensity of the b2 parameter, where for M4 b2 was lowest, producing a decelerating mortality instead of a typical logistic mortality. From the fully age-structured mortality rates, we calculated survival probabilities at every age interval [x, x + ∆x) as where ∆x = 1. We modelled age-specific fecundity rate b(x) using the following flexible exponential function with a quadratic effect as a function of age where γ0, γ1, γ2 > 0. In Eq. S3, γ is a vector of parameters such that γ0 controls the maximum offspring produced when x = γ2, γ1 determines how fast fecundity rate increases with age, and γ2 represents the age at which the maximum fecundity is achieved. As with mortality, we varied the values of γ to simulate a range fecundity trajectories. We calculated fecundity for a given age interval [x, x + ∆x) as With the resulting fecundity and survival probabilities, we constructed Leslie matrices (23) for these fully agestructured (Aa) trajectories and calculated the intrinsic population growth rate λ as the dominant eigenvalue together with the corresponding stable age distribution as the right eigenvector for λ, u = [u0, u1, . . . , uω] where we define ω as the age when only 0.1% of the population remains alive (i.e. S(x) = 0.001) (8).
We then constructed the corresponding two-stage models by assuming that juvenile survival, p0, p1, . . . , pα−1, was equivalent for both models with α representing the age at maturity, and calculated adult survival (pc) and fecundity (mc) for the two-stage model as the average survival and fecundity for adults (x ≥ α) in the fully age-structured model as respectively, where the deterministic transition matrices produce the same dominant eigen value, this is λa = λc. We calculated deterministic population growth rates, λa and λc (8), with function eigen() from the base package in the statistical software R (34).
B. Stochastic perturbations to demographic rates. To simulate stochastic perturbations for both matrices Aa and Ac, we simulated environmental shocks by randomly drawing values zt and vt from a normal distribution with mean 0 and variance 1. We then used a logistic function to perturb the age-specific survivals as where g(px) = log[px/(1 − px)], which limits px,t ∈ [0, 1]. Similarly, we used a logistic function for the shocks on fecundity as where h(bx) = log[bx/(M − bx)] and M = 2 is the maximum offspring, which limits the range of values for bx,t ∈ [0, M ]. We simulated three different scenarios: a) negative covariation between survival and fecundity (i.e. zt = −vt); b) no covariation between survival and fecundity (i.e. zt and vt independent); and c) positive covariation between survival and fecundity (i.e. zt = vt).
For the analysis of extinction times, we included demographic stochasticity following (15) as additional random shocks ht with mean 0 and variance inversely proportional to the population size, namely Var [ht] = σ 2 h /Nt, where we set σ 2 h = 1.
C. Comparing the distributions ofλa andλc. For both matrices, namely Aa and Ac, we produced stochastic perturbations such that, at time t the demographic rates are affected following the procedure we outlined in the previous section. We then constructed the corresponding stochastic matrices Aa,t and Ac,t. For simplicity here we refer to the stochastic matrices as At, noting that the same procedure was followed for the full age-dependent and the constant-adult matrices. Thus, at each iteration we estimate the population vector nt as where, at step t, we scale nt−1 to sum to unity as and thus Eq. S11 becomes wt = At wt−1. After this standardization, we calculate the population growth rate at time t as wx,t. [S13] We ran simulations of the perturbation analysis for T = 200 time steps, from which we calculated the averages We repeated these simulations 2 000 times for each combination of demographic rates and scenario. We then constructed the distributions ofλa andλc and calculated Kullback-Leibler discrepancies (KL) (21,26) between the resulting distributions. The KL discrepancies allowed us to determine the amount of information lost when estimating the stochastic population growth rate with a one-adult-stage model. To calculate the KL values we defined Pa(λ) as the distribution ofλa for a given combination of mortality and fecundity trajectories and under one of the three scenarios, and Pc(λ) the corresponding distribution forλc. The Kullback-Leibler discrepancy of Pc with respect to Pa, namely KL(Pa, Pc), is calculated as [S14] Thus, if both distributions are equal, then KL(Pa, Pc) = 0 and there is therefore no loss of information. As KL(Pa, Pc) increases, so does the loss of information.
To reduce the asymmetry between KL(Pa, Pc) and KL(Pc, Pa), and to make the results easier to interpret, we used McCulloch's (26) calibration. With this calibration, a KL(Pa, Pc) is equal to to a KL value between a Bernouilli distribution with probability (parameter) 0.5 and one with parameter q. If q is close to 0.5, this means that Pa and Pc are as close to each other as B(0.5) and B(q). The value of q is obtained as where k = KL(Pa, Pc). Thus, This calibration can be interpreted as a probability between 0.5 and 1, where 0.5 implies that both distributions are equal, and 1 that there is no overlap between them.

Review and alternative decomposition of E [λ t ] = λ e and Var [λ t ] = V λ
The following section provides alternative decompositions for the expected value and variance of λt. However, these are not intended as new results but as calculations that can help us explain the results from our simulations. These alternative decompositions are concurrent with the extensive previous work on the estimation of the expected value of λt and approximations to the long-run stochastic population growth rate (12,41,44,45), decompositions of the variance λt (6,7,13,42), and on the effect of the covariation between demographic rates on this variance (4,14,42,43).
A. Expressing λt in terms of demographic rates and age structure. We based our decompositions on the method proposed by Brown and Alexander (6) and Brown et al. (7), although as we mention above, we do not only focus on the variance of λt, but also on the expected value, which yields the expected population growth rate, λe. Furthermore, our method differs from theirs in that we decompose the variance in λt instead of the relative population growth rate δNi, while we do not take into account several within year mortality and recruitment episodes but assume that all mortality and recruitment can be summarized by a single yearly rate for each age. Under stochastic environments, the population growth rate is given by where Nt = ω x=0 nx,t is the total population size at time t and nx,t is the number of individuals in age x at time t, where ω is the maximum age the population can reach (42).
The equality in S16 can be reformulated as where bx ≥ 0 is the age-specific fecundity rate, px ∈ [0, 1] is the age-specific survival probability, while wx,t−1 = nx,t−1/Nt−1 is the proportions of individuals in age x at time t − 1 (for a full derivation see (42)). Thus, λt can be calculated as the sum of the proportion of individuals contributed by nx,t−1 to the next time t. We will note this proportion as yx,t = wx,t−1(bx,t + px,t). [S18] In addition, we define the expected population growth rate, noted as E [λt] = λe, as the theoretical mean of λt. This expected value can be approximated as for large t and where K is an initial burn-in period so that each simulated population reaches weak ergodicity in the age structure (11).
Our derivations below also require to define the expected values for the demographic rates as E [bx,t] = βx and E [px,t] = ρx for all x ∈ Ω, where Ω = {x ∈ N0 | x ≤ ω} is the set of discrete ages from 0 to the maximum age ω. These are the expected values of the age-specific fecundity rates and survival probabilities, respectively, while E [wx,t] = ηx is the expected values of the proportion of individuals in age x.
B. Expected value of λt. Here we provide a decomposition of λe based on the definition of λt in Eq. S17.
Let bx,t ≥ 0 and px,t ∈ [0, 1] be the age-and time-specific fecundity rate and survival probability for age x at time t of a given population, and let wx,t−1 ∈ [0, 1] be the proportion of individuals in age x at time t − 1, for x = 0, 1, 2, . . . , ω and t > 0, where ω is the maximum age the population can reach. If bx,t, px,t and wx,t−1, are age-and time-specific random variables with expected values given by E [bx,t] = βx, E [px,t] = ρx, and E [wx,t−1] = ηx, then the expected population growth rate, E [λt] = λe, will be Proof. Based on the results in Eqs. S17 and S18, we have that the expected value of λt can be expressed as the expectation of the sum of the elements in the vector of proportions y t = [y0,t, y1,t, . . . , yω,t], such that

Since by definition the covariance between any pair of random variables
However, under serially independent environments we have that C wb = Cwp = 0, which simplifies equation S20 to λe =λ. [S21] C. Variance of λt. Following the initial statement as in the previous section, whereby if bx,t, px,t and wx,t−1 are ageand time-dependent random variables for the fecundity rates, survival probabilities, and proportion of individuals in age x ∈ Ω, respectively, then the variance of the stochastic population growth rate, Var [λt] = V λ , is Therefore, to calculate the variance in λt it is necessary first to estimate the variance for the elements of the vector of proportions y t = [y0,t, y1,t, . . . , yω,t]. For simplicity of notation, we will not include the subindices t and t − 1. Thus, based on Eq. S18, we have that Expressing the variance in terms of its expectation and using the property Var [X] = E X 2 − E [X] 2 for any random varialbe X, we obtain Expanding the square we have that Next, we require providing the covariance between any two elements yi,t and yj,t of y, where i = j and i, j ∈ Ω. Once more, we will still omit the subindices t and t − 1. First, we can express the covariance between any pair of porportions yi and yj in terms of expectations as Then, we can replace yi,t and yj,t by their corresponding functions of demographic rates as in Eq. S18 which yields

Approximations to E [r t ] = r e
For any random variable X, with mean E [X] = µX and variance E (X − µx) 2 = σ 2 X , and a function f (X) n ≥ 2 times differentiable in the range of X, we have that the expectation E [f (X)] = re, known as the long-run population growth rate (45), can be approximated with the second order Taylor expansion given by Thus, for f (λt) = ln λt = rt, where E [λt] = λe calculated as in Eq. S20 and Var [λt] = V λ from Eq. S22, the second order Taylor approximation to E [rt] = re will take the form This formulation is concurrent with the small noise approximation to the long-run population growth rate proposed by Tuljapurkar (45) for serially uncorrelated environments, given by where λ0 is the dominant eigenvector of the matrix of average rates, and matrix τ 2 0 accounts for the variances of and covariances between the average demographic rates scaled by the sensitivities of λ0 to them. In other words, τ 2 0 provides an approximation to the variance in the population growth rate. The difference between Eqs. S27 and S28 is that the first is based on λe, which requires the estimation of the average demographic rates and average age-structure, while the second only depends on the matrix of average demographic rates. Both require averaging over the full environmental series. Thus, we would like to stress that the approximation in Eq. S27 requires long-term stochastic simulations and, therefore, cannot be based on a deterministic calculation, while Tuljapurkar and Haridas (44) provided a simplified method to calculate the small-noise approximation in Eq. S28.