Explaining mortality rate plateaus
-
Edited by Kenneth W. Wachter, University of California, Berkeley, CA, and approved October 24, 2001 (received for review May 8, 2001)
Abstract
We propose a stochastic model of aging to explain deviations from exponential growth in mortality rates commonly observed in empirical studies. Mortality rate plateaus are explained as a generic consequence of considering death in terms of first passage times for processes undergoing a random walk with drift. Simulations of populations with age-dependent distributions of viabilities agree with a wide array of experimental results. The influence of cohort size is well accounted for by the stochastic nature of the model.
Fundamental studies of the aging process have of late attracted the interest of researchers in a variety of disciplines, linking ideas and theories from biochemistry to mathematics (1–3). Much of this recent activity is due to the possibility that one of the supposedly fundamental tenets of aging, namely the exponential growth of mortality rates proposed by Gompertz (4), may fail to describe the behavior of observed populations adequately. More specifically, studies using populations or “cohorts” of Saccharomyces cerevisiae, Caenorhabditis elegans, Drosophila, and humans demonstrate that mortality rates tend to level off and even decrease at later stages of life (5–9). Attempts have been made to explain these plateaus via parabolic hazard functions (10), age-dependent demographics (11, 12), and phenomenological bifurcation models (8). In this paper we propose a simple model that incorporates heterogeneous dynamics to explain the generic plateau in mortality rates commonly observed in large cohorts of organisms.
Consider a population of N organisms with a distribution of
viabilities, v
i ≥ 0, where
v
i = 0 signifies death. The dynamics of an
individual viability will be modeled as follows:
where ɛ > 0 is a constant drift, σ > 0 is the standard
deviation of the fluctuations, and χi(t) is an
uncorrelated Gaussian random variable with zero mean and unit standard
deviation. The linear decline of viability is justified by the
observation of a linear decline of physiological functions noted by
Strehler and Mildvan (13) as well as similar results in more recent
surveys (14). The stochasticity in the system may be related to the
competition for resources, phenotypic differences, local environmental
changes, or even stochastic gene expression (15) but does not
necessarily depend on heterogeneity in the initial genotypic
distribution. The inclusion of stochasticity at the individual level
implies that Eq. 1 may be considered a changing frailty
model as opposed to a fixed frailty model. A fixed frailty model
preserves any initial heterogeneity in v throughout each
individual lifespan (16).
Biologically, this model states that each organism drifts toward death but with low probability may occasionally increase its long-term chances for survival. The probability of dying at time t is equivalent to the probability of first passage time P(t|v 0) of a random walk that begins at v = v 0 and reaches the origin, v = 0, at time t. The likelihood of death is controlled by the relative strength of drift and fluctuations. In this paper we explain the basic mechanism associated with first passage time problems and then proceed to show how such a model captures the essential features of late life mortality plateaus.
Theory of First Passage Time Problems
Eq. 1 may be better understood by considering the
limits of vanishing noise and then vanishing drift. When σ → 0, the
organisms move in lock step toward an inevitable death. The hazard
rate, μ(t), may be solved at t > 0 for
any initial distribution of viabilities, n
0(v),
where N
0 is the total number of initial
organisms, and D(t) is the number of organisms that die at
time t,
For a uniform distribution of initial viabilities, 0 <
v
i(0) < 1, the hazard rate reduces to
In the limit of slow drift the hazard rate grows exponentially for
small t and continues to grow until the system is left
desolate at t = 1/ɛ. At intermediate times the
hazard rate is not increasing exponentially as one might expect from a
Gompertz model. Regardless of its precise form, the monotonic increase
of mortality rates as evidenced in this simple example leads to the
natural question of what causes mortality rates to plateau in
populations of fruit flies, yeast, and other organisms.
The first step in answering this question is to consider the other limit of Eq. 1, namely ɛ → 0, when fluctuations dominate the dynamics. In this regime an individual viability v i(t) follows a random walk that ends when v ≤ 0. Qualitatively the removal of individuals with v ≤ 0 is tantamount to increasing the average viability of the remaining cohort. With time the average viability should increase, and therefore the hazard rate should decrease. It is important to note that as the population dies off, it will become more susceptible to fluctuations and may exhibit an intermittent rise in hazard rate near complete elimination. This caveat notwithstanding, we begin to see why the combination of these two effects, drift and fluctuation, might give rise to just the sort of behavior observed in large-scale mortality studies.
In order to simplify the analytical calculation of hazard rates we
rewrite Eq. 1 in the case of continuous time,
where W
i(t) is a stochastic Wiener process
that satisfies
dW
i(t′)dW
j(t) =
δijδ(t′ − t)dt. The difference between
Eqs. 1 and 5 is partly a matter of analytical
convenience¶; however, the
dynamics of aging may be related intimately to reproductive schedules,
which would imply the bulk of aging takes place at discrete intervals
(e.g. cell division in yeast) as opposed to taking place via a
continuous process of senescence. The qualitative behavior of both
models is essentially equivalent.
The closed form expression for the probability of crossing the origin
at t given an initial viability v
0,
i.e., the probability of dying between t and t +
dt, for systems behaving according to Eq. 5 is
This probability distribution is known as the inverse Gaussian
distribution (18–21). Researchers have considered its importance in
the context of a broad range of lifetime problems (22); for example, it
has proved useful in analyzing marriage durations (23), species
extinctions (24), and even the shelf life of food (25). Other
stochastic models of mortality have been proposed (26, 27), but to our
knowledge none has recognized the importance of the inverse Gaussian
distribution in the context of late-life mortality plateaus.
Given an initial distribution of viabilities,
n
0(v), the hazard rate in the continuous model
(Eq. 5) can be written as
where
is the density of individuals dying as a function of age. For a
completely homogeneous population, with
v
i(t) = v
0 and
N
0 individuals, D(t) =
N
0
P(t|v
0). For arbitrary heterogeneous
populations, D(t) may be calculated analytically or computed
numerically. In the case of homogeneous populations, e.g.
v
0 = 1 for all individuals at t
= 0, the hazard rate can be expressed analytically (28) as
where τ ≡ v
0/ɛ,
τr ≡
v
/σ2, and H(x)
is the complementary error function (29).
Although the expression for hazard rate is complicated, the behavior is controlled by two time scales: τ, corresponding to the mean lifetime, and τr, corresponding to the mean time for an organism to change its viability by v 0 via fluctuations alone. The mean lifetime does not depend on σ, the parameter responsible for changing frailty, although the variance of the lifetime (as well as higher moments) do depend on both ɛ and σ.
In Fig. 1 we show theoretical hazard rate curves when τ = 10 and τr = 2, 50, and 150. When τr ≫ τ, the system is dominated by drift, and hazard rates will increase rapidly after a brief delay at the initial stages before leveling off. When τr ≪ τ, the system is dominated by fluctuations, and hazard rates will actually decline with age after a brief increase. Finally, when τr > τ but still of the same order of magnitude, the system will exhibit an initial increase in hazard rate followed by a plateau. A range of behavior from essentially monotonic growth to plateaus to plateaus followed by declines is possible. Chhikara and Folks (30) have reached similar conclusions for the “failure rate” of any system, the probability of first passage times of which is modeled by Eq. 6. However, the rapid increase of hazard rate in all cases is slower than exponential. As it stands, the model is insufficient to explain mortality curves that have sustained Gompertz-type behavior at intermediate times. It is worth noting that any decrease in ɛ leads to an increase in the mean survival time, τ, as well as a postponement of the onset of the plateau, consistent with recent work on human populations (31).
Mortality rates according to Eq. 9 with τ = 10 and τr = 2, 50, and 150 corresponding to the solid, dashed, and dot-dashed curves, respectively.
One important point is that the asymptotic hazard rate for the inverse Gaussian distribution is nonzero, i.e., limt→∞ μ(t) > 0. However, for finite size populations the fluctuations inevitably lead to the eventual decline of the entire population, and therefore μ(t) = 1 at some finite value of t. The results of numerical simulations of Eq. 1 can be found in Fig. 2. Notice that the generic shape of the curves are the same as those found in Fig. 1, although the presence of fluctuations becomes important as the cohort dwindles in number.
Mortality rates obtained via numerical simulations of a homogeneous
population of N = 106 organisms, the
dynamics of which is that of Eq. 1, with
v
0 = 1, τ = 10, and
τr = 2, 50, and 150 corresponding to the solid,
dashed, and dot-dashed curves, respectively. Note that ɛ = 1/τ
and σ = τ
.
Because the presence of plateaus is related in some sense to the
difference between the most- and least-fit individuals in the cohort,
it is useful also to calculate the ratio between the viability of those
individuals closest and farthest from death. If we denote
v̄
as the average viability of the
top x fraction of individuals and
v̄
as the average viability of the
bottom x fraction of individuals, then the ratio of the most
fit to least fit can be expressed as r
x =
v̄
/v̄
. In the
stochastic simulations of Fig. 2, the ratio between the top 10% to the
bottom 10% leveled off at approximately r
0.1 =
15, 30, and 40 in the case of τ = 10 and
τr = 2, 50, and 150, respectively. Even if such
analysis does not constitute a formal proof, it does provide reassuring
evidence that the fitness gap between the most and least fit does not
grow without bound in this model.
Application to Mortality Rate Data
The usefulness of this simple model is to see whether experimental
data, specifically those exhibiting mortality plateaus, match with the
predictions of theory. Given a set of times X
i
specifying the lifespan of N organisms, the maximum
likelihood estimates (MLEs) of the parameters ɛ and σ in Eq.
6 are
and
where v
0 = 1 may be specified a
priori. Using these estimates of parameters we now compare
empirical observations to the analytical hazard rate assuming a
homogeneous distribution as well as to the results of simulations. The
empirical data on which we concentrate for the remainder of the paper
is from an experiment where N
0 = 1,203,646
genetically identical flies were observed for nearly 6 months, and the
number of deaths was recorded on a daily basis (5).
The MLE estimates of the parameters are ɛ̂ = 0.0448 and
σ̂ = 0.0975, which corresponds to τ̂ = 22.3 and
τ̂r = 105. In Fig.
3 we compare theoretical predictions for
the survival probability,
and the probability distribution P(t) in Eq.
6 to those calculated from the experimental data set. Notice
the good agreement not only with the survival probability (which
essentially is a cumulative measure) but with the actual probability
distribution of dying at a given age.
Survival probability at time t for experimental cohorts and theoretical predictions based on Eq. 6 where ɛ̂ = 0.0448 and σ̂ = 0.0975. (Inset) A comparison of the probability of dying at time t between experiment and theory. In both cases open circles signify data points from Carey et al. (5), and solid lines are predictions from theory.
By using these MLE estimates as a basis for further computational studies, we simulated Eq. 1 using v 0 = 1 and calculated the mortality rates for various size cohorts. In the study by Carey et al. (5), mortality rate is plotted against time for cohorts of different sizes. Not surprisingly the mortality rate is ill defined for small cohorts and well defined for large cohorts until mortality rate levels off and then exhibits strong fluctuations about a mean. To compare against experimental data we generated random cohorts by sampling from the fly mortality data of Carey et al. (5). Fig. 4 demonstrates how stochasticity and drift lead to both plateau- and size-dependent fluctuations very much in agreement with real data.
Age-specific mortality rates for cohorts of sizes ranging from n = 25 to 100,000. Going from top to bottom, rows 1 and 3 are randomly sampled from experimental data. Rows 2 and 4 are from simulations of Eq. 1 using MLE estimates as explained in Application to Mortality Rate Data.
The agreement should make it clear that small cohorts essentially veil the presence of mortality plateaus that would otherwise be evident in larger populations.
Discussion
This model of aging is in many ways an ideal construct. There is no upper limit on viability, no time-dependent drift to reflect developmental stages or the influence of natural selection, and no organism–organism interaction terms which might be important in large cohorts with limited resources. Modifications to the model perhaps could lead to better agreement with data but would not alter the main point we are trying to emphasize. The models we propose, Eqs. 1 and 5, and the behavior they exhibit demonstrate that mortality plateaus are a generic consequence of classifying lifespan as a first passage time problem for a random walk with drift. It is important also to note that although generic, mortality plateaus are not an inevitable outcome of this model for finite size cohorts.
There has been a long-standing interest on the effect of heterogeneity on mortality rates, although models have relied typically on notions of “fixed” frailty where individuals are fixed at the same relative risk compared to others in the same cohort (16). In the context of the present model, fixed frailty is analogous to having v 0 differ among individuals but change by drift alone. However, a review of gerontological studies shows that the variability of a majority of biological and cognitive health indicators tends to increase with age (32). Likewise a longitudinal study of blood pressure among nearly 4,000 men over a 40-year period found increasing variability with increasing age (33). These empirical results on age-related variability are more compatible with a changing frailty model even if some combination of changing and fixed frailty models is ultimately necessary to explain gerontological data.
It is only in the past half-century that substantive theories have been proposed to explain the underlying causes as well as the definition of aging. Different aging mechanisms are thought to involve oxidative damage from free radicals (34, 35), senescence genes (36), shortening telomeres (37), programmed cell death (38), simple “clogging” of cellular machinery in the nucleolus from hyperreplication of rDNA plasmids (39), mitotic misregulation (40), and of course combinations of all these as well (3, 41). The implications of these and other explanations of aging for population level mortality rates have not been well established. Subsuming these various mechanisms into a single measure of viability leads to the natural question of how ɛ and σ are correlated with body size, environmental conditions, and even genetic differences, if at all.
In order to test the model one might attempt a series of experiments wherein a deleterious treatment (e.g. radiation, DNA damaging agents, etc.) could be used to artificially control the level of environmental fluctuations influencing viability. For example, one might consider exposing a cohort of Drosophila to normally distributed levels of treatment with a nonzero level of treatment as control. In this way it would be possible to test the effects of increasing fluctuations σ with a fixed ɛ on cohort mortality rates.
Although this model demonstrates why an aging process should exhibit mortality plateaus, it does not claim to uncover the fundamental causes of aging. Bridging the gap between the causes and consequences of aging would provide a natural framework to explain how different mechanisms of aging lead to different types of mortality curves.
Acknowledgments
We thank L. Demetrius, P. Dodds, D. Rothman, N. Schorghofer, D. Sinclair, and K. Wachter for helpful comments. We also thank the anonymous member editor and referees for their suggestions and for pointing out a number of useful references. J.S.W. acknowledges the support of a National Defense Science and Engineering Graduate Fellowship. This work was supported in part by National Science Foundation Grant DMS-9709607.
Footnotes
-
↵ † To whom reprint requests should be addressed at: Department of Physics, Building 54-621, Massachusetts Institute of Technology, Cambridge, MA 02139. E-mail: jsweitz{at}segovia.mit.edu.
-
↵ § Present address: Department of Molecular and Cell Biology, University of California, Berkeley, CA 94709.
-
This paper was submitted directly (Track II) to the PNAS office.
-
↵ ¶ The probability of first passage times for the discrete version of the model are referred to as ruin problems (17).
- Abbreviation:
- MLE,
- maximum likelihood estimate
- Copyright © 2001, The National Academy of Sciences
