The long lives of primates and the ‘invariant rate of ageing’ hypothesis

Data for non-human primates

We obtained 30 datasets for six genera of non-human primates: sifaka (Propithecus spp), gracile capuchin monkey (Cebus spp), guenon (Cercopithecus spp), baboon (Papio spp), gorilla (Gorilla spp), and chimpanzee (Pan troglodytes) (Supplementary Data 1). Of these, 17 datasets correspond to long-term projects in the wild, while 13 were contributed by the non-profit Species360 from ZIMS¹⁸, which is the most extensive database of life history information for animals under human care.

Basic demographic functions

Let X be a random variable for ages at death, with observations x ≥ 0, and let μ (x|θ) be a continuous, non-negative parametric hazards rate or mortality function defined as

given that the limit exists, where ({boldsymbol{theta }}in {{mathbb{R}}}^{p}) is a p-dimensional vector of mortality parameters. The cumulative hazards rate is

which results in the survival function

The Cumulative distribution function (CDF) of ages at death is F (x | θ) = 1 – S (x | θ), and the probability density function (PDF) of ages at death is f (x | θ) = μ (x | θ) S (x | θ), for x ≥ 0. The remaining life expectancy after age x is calculated as

which yields a life expectancy at birth given by

The lifespan inequality at birth, as proposed by Demetrius^16,36 and later by Keyfitz¹⁷, is given by

Following Colchero et al.¹³, we define the lifespan equality as

For simplicity, henceforth we note the life expectancy, lifespan inequality and lifespan equality at birth as e(0 | θ) = e, H (0 | θ) = H, and ε (0 | θ) = ε, respectively.

Survival analysis

To estimate age-specific survival for all the wild populations of non-human primates, we modified the Bayesian model developed by Colchero et al.¹³ and Barthold et al.³⁷. This model is particularly appropriate for primate studies that follow individuals continuously within a study area and when individuals of one or both sexes can permanently leave the study area (out-migration), while other individuals can join the study population from other areas (in-migration). Thus, it allowed us to make inferences on age-specific survival (or mortality) and on the age at out-migration.

Here we use the five parameter Siler mortality function²⁵, as in Eq. (1) where θ = [a₀, a₀, c, b₀, b₁] is a vector of parameters to be estimated, and where a₀, b₀ ({mathbb{in }}{mathbb{R}}) and a₁, c, b₁ ≥ 0. For all species we studied, individuals of one or both sexes often leave their natal groups to join other neighbouring groups in a process commonly identified as natal dispersal. For some species, individuals who have undergone natal dispersal can then disperse additional times, described as secondary dispersal. Although dispersal within monitored groups (i.e. those belonging to the study area) does not affect the estimation of mortality, the fate of individuals that permanently leave the study area to join unmonitored groups can be mistaken for possible death. We identify this process as “out-migration”, which we classify as natal or immigrant out-migration, the first for natal and the second for secondary dispersals to unmonitored groups. This distinction is particularly relevant because not all out-migrations are identified as such, and therefore the fate of some individuals is unknown after their last detection. For these individuals we define a latent out-migration state at the time they were last detected, given by the random variable indicator O, with observations o_ij ∈ {0,1}, where o_ij = 1 if individual i out-migrated and o_ij = 0 otherwise, and where j = 1 denotes natal out-migration and j = 2 for immigrant out-migration. For known out-migrations, we automatically assign o_ij = 1. The model therefore estimates the Bernoulli probability of out-migration, π_j, such that O_ij ~ Bern(π_j). Those individuals assigned as exhibiting out-migration, as well as known emigrants and immigrants, contribute to the estimation of the distribution of ages at out-migration. Here, we define a gamma-distributed random variable V for ages at out-migration, with realisations v ≥ 0, where V_j | O_j = 1 ~ Gam(γ_j1, γ_j2) and where γ_j1, γ_j2 > 0 are parameters to be estimated with j defined as above. The probability density function for the gamma distribution is g_V(v | γ_j1, γ_j2) for v ≥ 0, with v = x_l – α_j, where x_l is the age at last detection and α_j is the minimum age at natal or immigrant out-migration.

In addition, since not all individuals have known birth dates, the model samples the unknown births b_i as x_il = t_il – b_i, where t_il is the time of last detection for individual i. The likelihood is then defined as

where x_if is the age at first detection, given by x_if = t_if – b_i, with t_if as the corresponding time of first detection. The parameter vectors γ₁ and γ₂ are for natal and immigrant out-migration, respectively. In other words, individuals with o_ij = 0 are assumed to have died shortly after the last detection, while those with o_ij = 1 are censored and contribute to the estimation of the distribution of ages at out-migration. The full Bayesian posterior is then given by

where the first term on the right-hand-side of Eq. (10) is the likelihood in Eq. (9), and the following terms are the priors for the unknown parameters. The vector π = [π₁, π₂] is the vector of probabilities of out-migration while the subscripts u and k refer to unknown and known, respectively.

Following Colchero et al.¹³, we used published data, expert information and an agent-based model to estimate the mortality and out-migration prior parameters for each population. We assumed a normal (or truncated normal distribution depending on the parameter’s support) for all the parameters. We used vague priors for the mortality and natal out-migration parameters (sd = 10), and informative priors for the immigrant out-migration parameters (sd = 0.5). We ran six MCMC parallel chains for 25 000 iterations each with a burn-in of 5000 iterations for each population, and assessed convergence using potential scale reduction factor³⁸.

For the zoo data we used a simplified version of the model described above, which omitted all parts that related to out-migration. In order to produce Supplementary Figs. 1 and 2, we used the same method as for the zoo data on the human life tables. To achieve this, we created an individual level dataset from the l_x column of each population, and then fitted the Siler model to this simulated data. It is important to note that the Siler model provides a close fit to the nonhuman primate data and to high-mortality human populations, although it does not provide the best fit to low-mortality human populations, in part due to the late life mortality plateau common among human populations³⁹ (Supplementary Fig. 6). It is therefore possible that the values of the mortality parameter b₁ we report in Supplementary Data 2 for the human populations are under-estimated. Nonetheless, and for the purposes of our analyses, the Siler fits to the human populations we considered here are reasonable (Supplementary Fig. 6) and we can therefore confidently state that the limitations of the Siler model do not affect the generality of our results.

Estimation of life expectancy and lifespan equality

Based on the results of the Bayesian inference models, we calculated life expectancy at birth as

where S (x) is the cumulative survival function as defined in Eq. (4) and where (hat{{boldsymbol{theta }}}) is the vector of mortality parameters calculated as the mean of the conditional posterior densities from the survival analysis described above. We calculated the lifespan inequality^17,36, H, as

from which we calculated lifespan equality, ε, as in Eq. (8). We calculated both measures for each of the study populations, and performed weighted least squares regressions for each genus, with weights given by the reciprocal of the standard error of the estimated life expectancies.

Sensitivities of life expectancy and lifespan equality to mortality parameters

As we mentioned above, for simplicity of notation, we will express all demographic functions by their variable notation (e.g. e = e (0 | θ), S = S (x | θ), etc.), while we will alternatively note first partial derivatives, for instance the derivative of e with respect to a given mortality parameter θ ∈ θ, as e_θ or ∂e / ∂θ.

Proposition: If ({S:}{{mathbb{R}}}_{ge 0}to left[{mathrm{0,1}}right]) is a continuous non-increasing parametric survival function with parameter vector ({boldsymbol{theta }}{boldsymbol{in }}{{mathbb{R}}}^{{boldsymbol{p}}}), with continuous differentiable cumulative hazards function ({U:}{{mathbb{R}}}_{ge 0}to {{mathbb{R}}}_{ge 0}), and with life expetancy at birth, lifespan inequality and lifespan equality as in Eqs. (4)-(6), respectively, then the sensitivity of life expectancy, e, to a given parameter θ ∈ θ is

while the sensitivity of lifespan equality to θ is

where

is the sensitivity of the survival function at age x to changes in parameter θ.

Proof. The sensitivity of lifespan equality to changes in θ is derived from

which, by Leibnitz’s rule, Eq. (16) becomes

The sensitivity of lifespan equality to changes in θ can be calculated as

By the quotient and Leibnitz’s rules, Eq. (18) can be modified as

The first term in Eq. (19) can be further decomposed by the product rule, while the second term can be modified following the equality for H in Eq. (7), which yields

By the chain rule, we have that (frac{partial U}{partial theta }=-frac{partial }{partial theta }{{log }},S=-frac{1}{S}frac{partial S}{partial theta }), which modifies Eq. (20) as

hence completing the proof. ∎

Changes in parameters along the genus lines

From the results in Eqs. (13) and (14), we calculated the vectors of change (gradient vectors) at any point (leftlangle {e}_{j},{varepsilon }_{j}rightrangle) of the life expectancy-lifespan equality landscape, as a function of each of the Siler mortality parameters (See Fig. 2A, B).

To quantify the amount of change of each parameter along the genus lines, we derived the sensitivities of a given mortality parameter θ to changes in life expectancy and lifespan equality, namely (frac{partial theta }{partial e}=frac{1}{{e}_{theta }}) for ({e}_{theta }, ne, 0,) and (frac{partial theta }{partial varepsilon }=frac{1}{{varepsilon }_{theta }}) for ({varepsilon }_{theta }, ne, 0). With these sensitivities we calculated the gradient vector

for any parameter at any point along the genus lines. Here we find a linear relationship between life expectancy and lifespan equality, given by

for i = 1, …, n_k, where n_k is the number of populations for genus k, and ({hat{varepsilon }}_{{ik}})is the fitted value of lifespan equality for population i in genus k, and β_0k and β_1k are linear regresssion parameters for genus k. To estimate the amount of change in parameter θ along the line for genus k, we can solve the path integral

where path C_k is determined by the linear model for genus k and (d{bf{r}}=leftlangle {de},dhat{varepsilon }rightrangle =leftlangle {de},{d; m}left(eright)rightrangle) is the rate of change in the velocity vector ({bf{r}}=leftlangle e,hat{varepsilon }rightrangle =leftlangle e,mleft(eright)rightrangle).

In order to compare results between the different mortality parameters in vector θ, we use the transformation g(θ) = log θ, which yields the following partial derivatives

and

Thus the gradient vector becomes

while the path integral in Eq. (24) is modified accordingly. In short, the path integral ({Theta }_{j}) provides a measure of the relative change in parameter θ along the genus line (Fig. 3). To allow comparisons between all genera, we scaled the values of each path integral by the length of each line.

Applications to the Siler mortality model

The Cumulative hazards for the Siler mortality model in Eq. (7) is given by

The sensitivities in Eqs. (13) and (14) require calculating S_θ for all θ ∈ θ. Treating S (x) as the function composition W (V), where W = exp(x) and V = – U, then S_θ is

where U_θ is the first derivative of U(x | θ) with respect to θ. For each of the Siler mortality parameters, we then have

All analyses were performed in the free open source programme R⁴⁰. The R functions we created for this project can be found in⁴¹.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com