The long lives of primates and the ‘invariant rate of ageing’ hypothesis
Data for non-human primatesWe 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 ZIMS18, which is the most extensive database of life history information for animals under human care.Basic demographic functionsLet 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$$mu left(x,|,{boldsymbol{theta }}right)=mathop{{rm{lim}}}limits_{Delta xto 0}frac{{{Pr }}(x < Xle x+Delta x|X > x)}{Delta x},$$
(2)
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$$Uleft(x|{boldsymbol{theta }}right)=int_{0}^{x}mu (t|{boldsymbol{theta }}){dt},$$
(3)
which results in the survival function$$S(x|{boldsymbol{theta }})={{exp }}[-U(x|{boldsymbol{theta }})].$$
(4)
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$$eleft(x{rm{|}}{boldsymbol{theta }}right) = frac{{int }_{x}^{{{infty }}}{tf}(t|{boldsymbol{theta }}){dt}}{Fleft({{infty }}right)-Fleft(xright)}\ =frac{{int }_{x}^{{{infty }}}S(t|{boldsymbol{theta }}){dt}}{Sleft(xright)},$$
(5)
which yields a life expectancy at birth given by$$eleft(0{rm{|}}{boldsymbol{theta }}right)={int }_{0}^{{{infty }}}S(x|{boldsymbol{theta }}){dx}.$$
(6)
The lifespan inequality at birth, as proposed by Demetrius16,36 and later by Keyfitz17, is given by$$H(0|{boldsymbol{theta }}) =-frac{{int }_{0}^{{{infty }}}S(x|{boldsymbol{theta }}){{log }}[S(x|{boldsymbol{theta }})]{dx}}{e(0|{boldsymbol{theta }})}\ = frac{{int }_{0}^{{{infty }}}S(x|{boldsymbol{theta }})U(x|{boldsymbol{theta }}){dx}}{e(0|{boldsymbol{theta }})}.$$
(7)
Following Colchero et al.13, we define the lifespan equality as$$varepsilon (x|{boldsymbol{theta }})=-{{log }}[H(x|{boldsymbol{theta }})].$$
(8)
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 analysisTo estimate age-specific survival for all the wild populations of non-human primates, we modified the Bayesian model developed by Colchero et al.13 and Barthold et al.37. 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 function25, as in Eq. (1) where θ = [a0, a0, c, b0, b1] is a vector of parameters to be estimated, and where a0, b0 ({mathbb{in }}{mathbb{R}}) and a1, c, b1 ≥ 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 oij ∈ {0,1}, where oij = 1 if individual i out-migrated and oij = 0 otherwise, and where j = 1 denotes natal out-migration and j = 2 for immigrant out-migration. For known out-migrations, we automatically assign oij = 1. The model therefore estimates the Bernoulli probability of out-migration, πj, such that Oij ~ 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 Vj | Oj = 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 gV(v | γj1, γj2) for v ≥ 0, with v = xl – αj, where xl 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 bi as xil = til – bi, where til is the time of last detection for individual i. The likelihood is then defined as$$p({x}_{{il}},{x}_{{if}},|,{boldsymbol{theta}},{boldsymbol{gamma}}_{1},{boldsymbol{gamma}}_{2},{pi }_{j},{o}_{ij})=left{begin{array}{cc}frac{fleft({x}_{il}right)}{Sleft({x}_{if}right)}({1-pi }_{j})hfill& {text{if}}; o_{{ij}}=0\ frac{Sleft({x}_{{il}}right)}{Sleft({x}_{{if}}right)}{pi }_{j}{g}_{V}({x}_{{il}}-{alpha }_{j})& {text{if}}; o_{{ij}}=1end{array}right.,$$
(9)
where xif is the age at first detection, given by xif = tif – bi, with tif as the corresponding time of first detection. The parameter vectors γ1 and γ2 are for natal and immigrant out-migration, respectively. In other words, individuals with oij = 0 are assumed to have died shortly after the last detection, while those with oij = 1 are censored and contribute to the estimation of the distribution of ages at out-migration. The full Bayesian posterior is then given by$$pleft({boldsymbol{theta }}{boldsymbol{,}}{{boldsymbol{gamma }}}_{1},{{boldsymbol{gamma }}}_{2},{boldsymbol{pi }},{{bf{b}}}_{u},{{bf{o}}}_{u},|,{{bf{b}}}_{k},{{bf{o}}}_{k},{{bf{t}}}_{f},{{bf{t}}}_{l}right) propto ; pleft({{bf{x}}}_{l},{{bf{x}}}_{f},|,{boldsymbol{theta }},{{boldsymbol{gamma }}}_{1},{{boldsymbol{gamma }}}_{2},{boldsymbol{pi }},{bf{d}}right)\ , times pleft({boldsymbol{theta }}right)pleft({{boldsymbol{gamma }}}_{1}right)pleft({{boldsymbol{gamma }}}_{2}right)pleft({boldsymbol{pi }}right),$$
(10)
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 π = [π1, π2] is the vector of probabilities of out-migration while the subscripts u and k refer to unknown and known, respectively.Following Colchero et al.13, 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 factor38.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 lx 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 populations39 (Supplementary Fig. 6). It is therefore possible that the values of the mortality parameter b1 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 equalityBased on the results of the Bayesian inference models, we calculated life expectancy at birth as$$e= int_{0}^{{infty}}Sleft(t| {hat{boldsymbol{theta }}}right){dt},$$
(11)
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 inequality17,36, H, as$$H=-frac{1}{e}int_{0}^{{{infty }}}Sleft(x{rm{|}}hat{{boldsymbol{theta }}}right){{log }}left[Sleft(x|hat{{boldsymbol{theta }}}right)right]{dx},$$
(12)
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 parametersAs 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$${e}_{theta }=frac{partial e}{partial theta }=int_{0}^{{{infty }}}{S}_{theta }{dx},$$
(13)
while the sensitivity of lifespan equality to θ is$${varepsilon }_{theta }=frac{partial varepsilon }{partial theta }=frac{{e}_{theta }left(1+{H}^{-1}right)-{H}^{-1}{int }_{0}^{{{infty }}}{S}_{theta }{Udx}}{e},$$
(14)
where$${S}_{theta }=frac{partial }{partial theta }S(x|{boldsymbol{theta }})$$
(15)
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$${e}_{theta }=frac{partial }{partial theta }int_{0}^{{{infty }}}{Sdx},$$
(16)
which, by Leibnitz’s rule, Eq. (16) becomes$${e}_{theta }=int_{0}^{{{infty }}}frac{partial S}{partial theta }{dx}=int_{0}^{{{infty }}}{S}_{theta }{dx}.$$
(17)
The sensitivity of lifespan equality to changes in θ can be calculated as$${varepsilon }_{theta } =frac{partial }{partial theta }left(-{{log }}, Hright)\ =-frac{partial }{partial theta }{{log }}, H\ =-frac{1}{H}frac{partial H}{partial theta }\ =-frac{1}{H}frac{partial }{partial theta }left(frac{{int }_{0}^{{{infty }}}{SUdx}}{e}right).$$
(18)
By the quotient and Leibnitz’s rules, Eq. (18) can be modified as$${varepsilon }_{theta } =-frac{1}{H{e}^{2}}left[frac{partial }{partial theta }left(int _{0}^{{{infty }}}{SUdx}right)e-left(int _{0}^{{{infty }}}{SUdx}right)frac{partial e}{partial theta }right]\ =-frac{1}{{He}}int _{0}^{{{infty }}}frac{partial }{partial theta }left({SU}right){dx}+frac{1}{{He}}frac{int _{0}^{{{infty }}}{SUdx}}{e}frac{partial e}{partial theta }.$$
(19)
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$${varepsilon }_{theta } =-frac{1}{{He}}int _{0}^{{{infty }}}left(frac{partial S}{partial theta }U+Sfrac{partial U}{partial theta }right){dx}+frac{1}{e}{e}_{theta }\ =-frac{1}{{He}}left(int _{0}^{{{infty }}}{S}_{theta }{Udx}+int _{0}^{{{infty }}}Sfrac{partial U}{partial theta }dxright)+frac{1}{e}{e}_{theta }.$$
(20)
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$${varepsilon }_{theta } = , -frac{1}{{He}}left(int _{0}^{infty }{S}_{theta }{Udx}-int _{0}^{infty }frac{partial S}{partial theta }{dx}right)+frac{1}{e}{e}_{theta }\ = , -frac{1}{{He}}left(int _{0}^{infty }{S}_{theta }{Udx}-{e}_{theta }right)+frac{1}{e}{e}_{theta }\ = , -frac{int _{0}^{infty }{S}_{theta }{Udx}}{{He}}+frac{{e}_{theta }}{e}left(1+frac{1}{H}right)\ =, frac{{e}_{theta }left(1+{H}^{-1}right)-{H}^{-1}int _{0}^{infty }{S}_{theta }{Udx}}{e},$$
(21)
hence completing the proof. ∎Changes in parameters along the genus linesFrom 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$$nabla theta =leftlangle frac{partial theta }{partial e},frac{partial theta }{partial varepsilon }rightrangle$$
(22)
for any parameter at any point along the genus lines. Here we find a linear relationship between life expectancy and lifespan equality, given by$$mleft({e}_{{ik}}right)={hat{varepsilon }}_{{ik}}={beta }_{0k}+{beta }_{1k}{e}_{{ik}},$$
(23)
for i = 1, …, nk, where nk 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$${Theta }_{k}=int _{{C}_{k}}nabla theta d{bf{r}},$$
(24)
where path Ck 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$$frac{partial }{partial e}gleft(theta right)=frac{1}{theta }frac{partial theta }{partial e}$$
(25)
and$$frac{partial }{partial varepsilon }gleft(theta right)=frac{1}{theta }frac{partial theta }{partial varepsilon }.$$
(26)
Thus the gradient vector becomes$$nabla theta =leftlangle frac{partial }{partial e}gleft(theta right),frac{partial }{partial varepsilon }gleft(theta right)rightrangle$$
(27)
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 modelThe Cumulative hazards for the Siler mortality model in Eq. (7) is given by$$Uleft(xright)=frac{{e}^{{a}_{0}}}{{a}_{1}}left(1-{e}^{{-a}_{1}x}right)+{cx}+frac{{e}^{{b}_{0}}}{{b}_{1}}left({e}^{{b}_{1}x}-1right),$$
(28)
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$${S}_{theta }=frac{{dW}}{{dV}}{V}_{theta }=-S{U}_{theta },$$
(29)
where Uθ is the first derivative of U(x | θ) with respect to θ. For each of the Siler mortality parameters, we then have$${S}_{{a}_{0}}=S(x|{boldsymbol{theta }})frac{{e}^{{a}_{0}}}{{a}_{1}}left({e}^{-{a}_{1}x}-1right)$$
(30)
$${S}_{{a}_{1}}=S(x|{boldsymbol{theta }})frac{{e}^{{a}_{0}}}{{a}_{1}}left[frac{1}{{a}_{1}}-{e}^{-{a}_{1}x}left(x+frac{1}{{a}_{1}}right)right]$$
(31)
$${S}_{c}=-S(x|{boldsymbol{theta }})x$$
(32)
$${S}_{{b}_{0}}=S(x|{boldsymbol{theta }})frac{{e}^{{b}_{0}}}{{b}_{1}}left(1-{e}^{{b}_{1}x}right)$$
(33)
$${S}_{{b}_{1}}=S(x|{boldsymbol{theta }})left[{e}^{{b}_{1}x}left(frac{1}{{b}_{1}}-xright)-frac{1}{{b}_{1}}right].$$
(34)
All analyses were performed in the free open source programme R40. The R functions we created for this project can be found in41.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More
