Original Research ARTICLE
Optimal Versus Realized Trajectories of Physiological Dysregulation in Aging and Their Relation to Sex-Specific Mortality Risk
- 1Biodemography of Aging Research Unit (BARU), Social Science Research Institute, Duke University, Durham, NC, USA
- 2Groupe de Recherche PRIMUS, Department of Family Medicine, CHUS-Fleurimont, University of Sherbrooke, Sherbrooke, QC, Canada
- 3The Danish Aging Research Center, University of Southern Denmark, Odense, Denmark
- 4Department of Clinical Genetics, Odense University Hospital, Odense, Denmark
- 5Department of Clinical Biochemistry and Pharmacology, Odense University Hospital, Odense, Denmark
While longitudinal changes in biomarker levels and their impact on health have been characterized for individual markers, little is known about how overall marker profiles may change during aging and affect mortality risk. We implemented the recently developed measure of physiological dysregulation based on the statistical distance of biomarker profiles in the framework of the stochastic process model of aging, using data on blood pressure, heart rate, cholesterol, glucose, hematocrit, body mass index, and mortality in the Framingham original cohort. This allowed us to evaluate how physiological dysregulation is related to different aging-related characteristics such as decline in stress resistance and adaptive capacity (which typically are not observed in the data and thus can be analyzed only indirectly), and, ultimately, to estimate how such dynamic relationships increase mortality risk with age. We found that physiological dysregulation increases with age; that increased dysregulation is associated with increased mortality, and increasingly so with age; and that, in most but not all cases, there is a decreasing ability to return quickly to baseline physiological state with age. We also revealed substantial sex differences in these processes, with women becoming dysregulated more quickly but with men showing a much greater sensitivity to dysregulation in terms of mortality risk.
The aging process involves many physiological changes affecting the homeostatic state of the organism (1). Studies of biomarkers of aging thus have the potential to simultaneously shed light on the underlying aging process and to provide measurement tools for understanding how aging proceeds within an individual. While a large number of studies have examined changes in biomarkers during the aging process, almost all of these have considered markers one at a time and without reference to how optimal levels of the markers change with age [see, e.g., Ref. (2)]. However, it is clear on the one hand that optimal levels can and do change with age (3, 4), and on the other that information contained in biomarkers is highly redundant, non-linear, and complex (5–7). While some studies have succeeded in incorporating one or the other of these complexities, none have yet attempted both simultaneously. The need for new methodological developments that would analytically integrate the number of biomarkers and complex interrelationships among them is recognized in the literature (2).
Recent advances in biodemographic models of aging resulted in development of mathematical models of aging that are based on biological theory, incorporate several essential mechanisms of aging-related changes in organisms, and work with relevant socio-demographic and other information collected in longitudinal data on aging. These models, known as the stochastic process models (SPM) or the quadratic hazard models (4, 8–10), allow us to evaluate “hidden components” of aging-related changes including adaptive capacity, resistance to stresses, physiological norms, and effects of allostatic adaptation – all variables that play important roles in the processes of aging (11–15); their inclusion in the model is therefore important to understand regulatory mechanisms driving observed aging-related changes in physiological variables and their influence on risks of death or developing a chronic disease, as well as for evaluating a genetic component in such processes. The advantage of the SPM-based approaches is that “hidden components of aging” can be estimated from the data although corresponding variables are not directly observed in the study.
Recently, Cohen et al. (5) presented a novel approach for measuring physiological dysregulation via the joint distribution of multiple biomarkers based on calculation of a multivariate distance called the Mahalanobis distance (which we denote DM) (16, 17). Such a measure represents the deviation of a current physiological state of an organism from the “normal” physiological state. In contrast to previous approaches, such as those used to calculate allostatic load, DM does not depend on a choice of biomarkers known to be associated with poor health, and thus avoids a potential for circularity (18). The interpretation of DM as physiological dysregulation has been validated by showing that it (a) increases with age within individuals, (b) predicts mortality controlling for age, (c) gains predictive power as more variables are included, and (d) is not particularly sensitive to which markers are included, nor to their individual correlations with age (5). It has even been shown to function in a wild bird species (19). While individual trajectories of DM with age have been examined in relation to health outcomes (20), the modeling approach used lacked the advantages of SPMs listed above, notably the ability to quantify adaptive capacity and optimized versus realized trajectories. Differences between men and women have also not been characterized.
The primary goal of this paper is to understand individual trajectories of dysregulation during aging in terms of adaptive capacity, mortality risk, and male–female differences. To do this, we will implement the measure of multivariate distance (DM) in the SPM framework. Joining the two approaches (SPM, DM) in the same model can further increase our capacity to uncover important aging-related processes. First, DM can be used as a quantitative measurement of how an individual departs from the physiological norm and is thus directly informative about parameters to be estimated by SPM. Second, although successfully applied to different data in different settings in one- or two-dimensional cases (3, 4, 21–27), the SPM may face computational difficulties in applications to multidimensional data because of a large number of parameters to be estimated and a computationally intensive likelihood maximization procedure. DM is useful in application of the SPM because it allows us to work with multiple physiological variables in a one-dimensional model while still allowing us to estimate and interpret all components of this one-dimensional SPM in the same way as in the original model. The SPM originates from the Woodbury–Manton model (28) and complexities associated with multidimensional versions of such models are recognized, see, e.g., Ref. (29), so the current approach is an important step forward in this research area. This paper is the first one addressing the combined effect of multiple physiological variables on mortality in the SPM (8) framework. We carry this out by applying the model to the Framingham Original Cohort (FHS) data. We estimate the model separately for female and male FHS participants to investigate sex differences in the dynamic behavior of the measure of physiological dysregulation and its possible relationships to the observed sex-specific differences in mortality risks.
Data and Methods
Measure of Statistical Distance and Physiological Dysregulation
For a given set of physiological variables (biomarkers) represented by a column vector X measured in an individual at age t, X(t), DM is defined as (5):
where is a vector of means, S is the variance–covariance matrix for respective variables calculated from some population defining the “normal” state (which can be the same population or some other “reference” population), and T denotes transposition.
Mathematical Model of Age Dynamics of Physiological Variables, Aging, and Mortality
Recent advances in mathematical modeling of aging allow us to link the age dynamics of physiological variables and mortality risks and evaluate “hidden components” of aging-related changes incorporating several major concepts of aging in the structure of the model. These models, originating from the Woodbury–Manton’s random walk model (28), became known as the SPM (8), see Section “Introduction.” Theoretical background of the models with survival functions induced by stochastic covariates is presented in Ref. (28, 30, 31). These results are directly applicable to the SPM version by Yashin et al. (8). The random evolution of covariates with age or time (which is a natural assumption, for example, for the evolution of biomarkers affecting survival) is not explicitly modeled in standard approaches, such as logistic or Cox models, and it is captured in such stochastic models. This is a fundamental feature of such models, which is an appealing property for practical applications, but the presence of stochastic covariates requires computations of marginal hazards and survival functions from the conditional hazards. When the conditional hazard is a quadratic function of such random covariates [as in Ref. (8)], such expressions can be computed as shown in Ref. (32, 33). Below, we briefly describe the SPM (8) in which we implement the measure of statistical distance (Eq. 1).
Briefly, the SPM consists of two equations representing the dynamics of chosen variables as a function of age t and conditional mortality rate at age t given the vector of the variables measured at respective ages. Let Y(t) be a K-dimensional stochastic process representing a (K × 1) vector of variables at age t. The conditional hazard of death is specified as:
Here, μ0(t) is the “residual” or “baseline” hazard that represents the mortality risk, which would remain if the variables in Y(t) followed the “optimal trajectory” represented by a (K × 1) vector function f0(t). This baseline hazard μ0(t) models the effect of other risk factors not included in Y(t) that impact the risk of death. The age trajectory of variables f0(t) that minimizes the risk of death is referred to as the “physiological norm.” The matrix Q(t) is a (K × K) positive-definite symmetric matrix.
The dynamics of the stochastic process Y(t) is given by the following equation:
with initial condition Y(t0). Here, W(t) is a vector Wiener process with independent components [W(t) is assumed to be independent of the initial vector Y(t0)], which defines random paths of the variables in Y(t). A matrix of diffusion coefficients b(t) regulates how this “randomness” propagates to variability of trajectories of Y(t).
The vector function f1(t) (referred to as the “mean allostatic state”) describes the effect of allostatic adaptation (13). This state characterizes the allostatically modified set point for homeostatic regulation, that is, the physiological state that organisms are forced to maintain by the process of homeostatic adaptation, which represents averaged effects of complex interplay among factors controlled by the ontogenetic program, senescence, and long-term stresses exceeding the limits of the organism’s homeostatic regulation. This dynamic behavior [i.e., that Y(t) adapts to changes in the function f1(t)] in the model is possible due to the presence of the negative feedback mechanism in Eq. 3 with coefficients of homeostatic regulation given by a matrix a(t). Age-related changes in these coefficients characterize changes in adaptive capacity with age. Specifically, the elements of this matrix characterize the rate of the adaptive response for any deviation of the variables from the state f1(t). This allows us to test hypotheses on decline in adaptive capacity with age [see also more discussion and examples in Ref. (8, 10)]. Note that the general model with age-dependent components in Eqs 2 and 3 can be simplified to the version with age-independent components. This may be helpful in some applications where testing the hypotheses on age-dependence of respective characteristics is not important but in our applications age-dependence of all components [except f0(t) and b(t)] is essential as it allows for addressing the hypotheses formulated below in Section “Application to Framingham Heart Study Data.”
The SPM described above can be applied to any set of variables for which longitudinal measurements are available in the data. The measure of physiological dysregulation (DM) represents an example of such a variable: it can be calculated for each individual at each observation point (exam) so that each individual has a trajectory of measurements of DM at different ages.
Mathematical Model of Age Trajectories of Physiological Dysregulation, Aging, and Mortality
Let Y(t) = DM[X(t)] as given by Eq. 1. This is now a one-dimensional process represented by Eq. 3 where respective components [W(t), f1(t), a(t), b(t)] are scalar counterparts of corresponding vectors/matrices. The one-dimensional version of Eq. 2 is:
where f0(t) and Q(t) are also scalars. As DM measures physiological dysregulation in an organism represented as deviations from the “normal state” (as specified in its calculation), it is natural to assume that f0(t) = 0 in Eq. 4. That is, in this implementation we assume that the “normal” values of physiological variables from the definition of DM (Eq. 1) minimize the risk of death at all ages [see also Section “Results and Discussion” regarding the specification of the “normal state” in DM calculations and alternative implementations of the “norm” f0(t)]. Note that the definition of DM involves calculations of a vector of means and the variance–covariance matrix for respective variables from some “reference” population defining the “normal” state. The question of which “reference” population should be used for these purposes can be dictated by data availability and research paradigms in specific studies, but the results are generally not overly sensitive to specifics of the choice (34). Reference populations that differ markedly from study populations across many demographic aspects (e.g., age, race, and sex) perform poorly, but differences in any single aspects have mostly minor effects. In our application, we used the population of individuals from the original Framingham cohort aged 40 years and younger at the baseline exam, as described below in Section “Application to Framingham Heart Study Data.”
Application to Framingham Heart Study Data
The Framingham Original Cohort was launched in 1948 and has continued with biennial examinations to the present. The FHS Original Cohort included 5,209 respondents (nearly all subjects were Caucasians) aged 28–62 years at baseline residing in Framingham, MA, USA, between 1948 and 1951. Data on 5,079 individuals (2,785 females, 2,294 males) from the Original Cohort were available for this study. Data relevant to our analyses include ages at exams, sex, ages at death/censoring (available data contain information on the number of days since the first exam until the event/censoring from which we calculated respective ages dividing the number of days by 365.25), and various physiological variables measured at exams. For our analyses, we selected physiological variables measured at a sufficient number of exams and whose dynamic characteristics have been shown to be related to the risk of death or onset of aging-related diseases in participants of the FHS Original Cohort in our earlier studies including those which used the original SPM approach [see, e.g., Ref. (3, 4, 24, 26, 27, 35)]. The list includes blood glucose (BG), body mass index (BMI), total cholesterol (CH), diastolic blood pressure (DBP), hematocrit (HC), systolic blood pressure (SBP), pulse pressure (PP), and ventricular rate (VR). Their relevance to research on aging has been discussed in the literature, see, e.g., an earlier publication by Manton et al. (36).
The original data on physiological measurements contain missing values (either by design when at some exam the respective variable was not measured or intermittent missing values when some individuals missed an exam). Following Engels and Diehr (37), we imputed intermittent missing longitudinal values for an individual using available data for that individual. Intermittent missing values were imputed using a linear approximation of available observations in the adjacent exams. For missing data at initial exams, two methods were used: the “next observation carried backward” method and the average of measurements at the next two exams. We repeated all computations using the data generated by these two methods and the results were not sensitive to the imputation method so only the results that used the latter method are reported here.
The data on original variables (BG, BMI, CH, DBP, HC, PP, SBP, and VR) were each transformed using the Box–Cox transformation and then standardized to a zero mean and a unit variance so that all transformed variables are on the same scale in the analyses. We used these standardized transformed variables in calculations of DM, as defined in Eq. 1. We calculated the means and the variance–covariance matrix in Eq. 1 using baseline measurements for individuals aged 40 years and younger (altogether, there were 2,012 such individuals in the Original Cohort, 1,104 females and 908 males). All analyses were performed separately for females and males (to avoid any assumptions on relationships between parameters of the model in two sexes) using the sets of physiological variables included in calculations of DM variants, as described below.
The values of DM depend on the specific set of variables included in its definition. Therefore, to test if the results are sensitive to the choice of variables in DM, it is necessary to perform analyses with different variants of DM based on different subsets of variables. The most comprehensive approach would be to use all possible combinations of variables as in Cohen et al. (5). However, as the present study applies the SPM that involves intensive computations, the analyses of all such variants are not feasible. Therefore, we selected a “basic” variant of DM and performed sensitivity analyses for an additional (limited) set of variants to test whether the results are sensitive to the choice of the variables in the construction of DM. Our “basic” variant is “the most complete” DM when all variables are included in its construction. Note that, as PP = SBP − DBP, we cannot use these three variables simultaneously in analyses. Therefore, we considered three “basic” variants each including seven variables: BG, BMI, CH, HC, VR, and two variables out of DBP, PP, and SBP. We also computed several additional variants of DM to perform sensitivity analyses to test: (a) whether the removal of one variable from such “basic” variants substantially changes the results; and (b) whether the removal of HC and VR (which, unlike the other variables, were not measured until exam 4) from such “basic” variants substantially changes the results. Altogether, we analyzed 24 different DM variants.
In our applications of SPM (Eqs 3 and 4) to these DM variants we used the following functional forms of the model’s components: (1) a linear function of age for the feedback coefficient: a(t) = aY + bY(t − tmin), where aY < 0, bY ≥ 0, and tmin = 28 (which is the minimal age at the first exam in the FHS Original Cohort); (2) a linear function of age for the mean allostatic trajectory: ; (3) a linear function of age for the multiplier in the quadratic hazard term [named “vulnerability index” in Arbeev et al. (4) as it characterizes the “robustness,” or “vulnerability,” a component of stress resistance]: Q(t) = aQ + bQt; (4) the Gompertz baseline hazard: ; (5) a constant diffusion coefficient: b(t) = σ1; and (6) a normally distributed initial value .
This is a parsimonious specification of the model that still allows us to test several hypotheses on the impact of systemic dysregulation on mortality risk and its relation to different aging-related processes represented by components of the model:
(1) H0: Q(t) = 0 [i.e., aQ = 0 and bQ = 0, so that there is no quadratic term in the hazard rate (Eq. 4) and mortality is described by the baseline Gompertz rate μ0(t)]. This is the most important hypothesis because the failure to reject this hypothesis would not allow us to claim that DM is a legitimate characteristic for capturing physiological dysregulation affecting mortality risk in this sample.
(2) H0: Q(t) = aQ (i.e., bQ = 0). Rejection of this hypothesis would mean that the effect of physiological dysregulation on mortality risk is age dependent. This hypothesis allows us to make inferences on aging-related change in stress resistance associated with the variables used in DM. For example, if Q(t) increases with age then the J-shape of the hazard rate (as a function of DM at a fixed age t) narrows with age indicating that the same value of DM induces a larger increase in mortality risk at old ages compared to younger ages, which can be interpreted as the aging-related decline in resistance to stresses associated with the dynamics of the respective variables, see Yashin et al. (8) and Arbeev et al. (4).
(3) H0: f1(t) = 0 (i.e., and ). This hypothesis allows us to make inferences on the presence of a systemic dysregulation for the specific set of variables used in DM. Rejection of this hypothesis would mean that there is a systemic dysregulation in an organism that forces the trajectories of physiological variables to deviate from their “normal values” specified in calculation of DM [or, equivalently, forces Y(t) to deviate from f0(t)].
(4) H0: (i.e., ). This hypothesis allows us to make inferences on the changes in the level of systemic dysregulation with age for the specific set of variables used in DM. Rejection of this hypothesis would mean that the level which the trajectory of Y(t) is forced to follow changes with age. For example, increasing f1(t) with age would mean an increasing level of systemic physiological dysregulation with age (so that the trajectories of physiological variables deviate further with age from their “norms” defined at younger ages), which can be associated with the manifestation of the aging process.
(5) H0: a(t) = aY (i.e., bY = 0). This hypothesis allows us to make inferences on the decline in adaptive capacity with age for the specific set of variables used in DM. Rejection of this hypothesis would mean that the value of the feedback coefficient in Eq. 3 changes with age. As noted, this feedback coefficient in the model is associated with the adaptive capacity of an organism [i.e., the rate of the adaptive response for any deviation of Y(t) from f1(t)]. In general, the larger the absolute value of a(t), the faster Y(t) tends to f1(t). This means that if the absolute value of a(t) declines with age then more time is needed for the trajectory of Y(t) to go back to f1(t) at old ages compared to younger ages, which is the manifestation of the decline in adaptive capacity with age. Note also that Y(t) can deviate from f1(t) in both directions, i.e., closer to zero as well as toward larger values, and if the absolute value of a(t) declines with age then such periods when the trajectory of DM departs farther from f1(t) tend to be longer at old ages. This may mean that aging-related changes observed in the aging human body may also manifest effects of compensatory adaptation which tend to reduce dysregulation effects when f1(t) becomes large enough. These considerations stimulate further development of methods of dynamic modeling of aging-related changes and their connections with health and survival outcomes. See also further discussion on relationships between adaptive capacity, dysregulation and different components of aging-related changes in Section “Results and Discussion.”
Testing the above hypotheses involves fitting the restricted models (with restrictions on the parameters specified as in the parentheses above) along with the original (unrestricted) model using the same likelihood estimation procedure (8). As we deal with nested models in all cases, the likelihood ratio test can be used to make statistical inference (see Table S1 in Supplementary Material). The likelihood maximization in the SPM is performed using the constrained optimization procedure (implemented, e.g., in the MATLAB’s Optimization toolbox). Constrained optimization is needed because constraints on parameter values are necessary in respective components of the model both for the mortality risk as well as for the risk factor Y(t). The details can be found elsewhere [e.g., Ref. (8, 9)]. Also note that Y(t) can always be transformed (e.g., using the logarithm or the Box–Cox transformation) similar to the original variables X(t), if necessary. The original scale was used in the paper for the sake of interpretability of the resulting trajectories.
Results and Discussion
Estimates of Parameters and Results of Testing Hypotheses
Estimates of parameters of the SPM applied to different variants of DM calculated for individuals from the Framingham Original Cohort as described in Section “Application to Framingham Heart Study Data” are given in Table 1 (for females) and Table 2 (for males). The tables show that, although the parameter estimates differ for different variants of DM, they still follow a common pattern so that generally the components of the model look similar in applications to considered DM variants (a few exceptions are described below). The tables also reveal systematic differences between estimates of parameters in females and males. This is addressed in more detail in Section “Sex-Specific Differences in Estimates in the Model.”
Table 1. Estimates of parameters of the stochastic process model applied to different variants of DM calculated for females from the Framingham Original Cohort.
Table 2. Estimates of parameters of the stochastic process model applied to different variants of DM calculated for males from the Framingham Original Cohort.
Tables 1 and 2 also present the results of the hypothesis tests specified in Section “Application to Framingham Heart Study Data” (see notes in the tables for definitions of symbols or the absence of those used in different columns of the tables). Most importantly, the tables show that the null hypotheses H0: Q(t) = 0 is rejected for all variants of DM (p < 0.0001 for all DM in males and most DM in females, with maximal p = 0.0007 for DM with BG, BMI, CH, HC, SBP, and VR). This allows us to conclude that there is a quadratic term in the hazard rate (Eq. 4) so that mortality is not entirely captured by the baseline Gompertz rate μ0(t), regardless of the variant of DM. This means that DM captures effects of deviations of physiological variables on mortality in this application. All considered DM affect mortality risk in this sample so that non-zero values of DM result in a higher mortality rate compared to the baseline mortality rates at respective ages. Note also that, although we have non-zero estimates of the multiplier Q(t) and the quadratic term is present in Eq. 4, we also have non-zero estimates of the baseline mortality rate μ0(t). There are many more risk factors that affect mortality risk and which are not included in definitions of DM and, correspondingly, in the quadratic part of the hazard. This “residual” mortality rate μ0(t) summarizes the effects of such unspecified factors affecting mortality risk.
The null hypothesis H0: Q(t) = aQ is also rejected in all cases, with p ≤ 0.0001 for males and p-values ranging between 0.0005 and 0.01 for females (except DM with BG, BMI, CH, DBP, and SBP for which p = 0.038). Tables 1 and 2 also show that in all cases the parameter bQ is positive. These observations indicate that the effect of physiological dysregulation (represented by DM) on mortality risk is age dependent and that the J-shape of the mortality rate (as a function of DM at any fixed age) narrows with age. This narrowing pattern with age has the effect that the same level of physiological dysregulation (i.e., the same value of DM) induces a larger increase in mortality risk at old ages than it does at younger ages (compared to the baseline mortality rate for respective age). This can be interpreted as the aging-related decline in resistance to stresses associated with deviant dynamics of respective physiological variables [see also Yashin et al. (8) and Arbeev et al. (4)].
The null hypotheses H0: f1(t) = 0 and H0: are rejected in all cases (all p < 0.0001). This result shows the presence of a systemic dysregulation in an organism that forces the trajectories of physiological variables to deviate from their “normal values” specified in calculation of DM (equivalently, forces DM to deviate from zero). The effect of this dysregulation is that, on average, an individual tends to have the values of physiological variables that deviate from the “normal” state (i.e., he/she has non-zero DM) so that the resulting mortality risk is higher than it could be if that person managed to keep his/her physiological variables equal to the “norm” (i.e., to have zero DM) in which case the mortality risk would be equal to the baseline mortality at the respective age. Note also that in applications to all variants of DM parameter is positive. This indicates an increasing level of systemic physiological dysregulation with age that manifests itself in trajectories of physiological variables drifting further away from their “norms” as age progresses.
The results for the null hypothesis H0: a(t) = aY are mixed. For many variants of DM (16 for females and 6 for males), the decline in the absolute value of the feedback coefficient a(t) with age (associated with the decline in adaptive capacity of an organism) is significant with p < 0.0001, and for some other variants, the decline is not significant, with some estimates bY = 0. This means that the conclusions about the decline in adaptive capacity depend on the variables used in calculations of DM and that different physiological variables have different behaviors in terms of the change with age in the ability of an organism to push them back to the trajectories specified by the mean allostatic trajectories. For those variables that indicated the decline in the absolute value of the feedback coefficient a(t) with age, this ability worsens with age and more time is needed for the trajectories of physiological variables to go back to their mean allostatic trajectories at old ages compared to the time needed at younger ages. However, for some physiological variables, this ability of an organism seems not to be affected by age. These results confirm our earlier observations on mixed patterns of age dynamics of adaptive capacity in analyses of different physiological variables and different outcomes (4, 26). A strong negative correlation between intercept and slope for the feedback coefficient a(t) (−0.86 for males, −0.92 for females) indicates that for the DM variants with worse initial adaptive capacity (i.e., smaller absolute values of aY) the subsequent decline tends to be slower than that for the DM variants with better initial adaptive capacity (i.e., larger absolute values of aY).
Sex-Specific Differences in Estimates in the Model
Tables 1 and 2 reveal systematic differences between estimates of parameters in females and males. These differences can be better understood from Figures 1 and 2. Figure 1 compares the estimates of female and male patterns of different components of the model applied to different variants of DM. It shows that the baseline mortality rate μ0(t) is higher in males. Males also have higher levels of the multiplier in the quadratic hazard term Q(t), which also increases faster with age than that in females. This means that males have a narrower J-shape of the mortality risk as a function of DM at each age and, moreover, this J-shape narrows even further at a faster rate in males than the J-shape for females. This confirms our earlier observations (4) and suggests that males have generally lower resistance to stresses associated with deviant dynamics of the respective physiological variables than females: at each age males are more vulnerable to deviations from the “normal” state (i.e., the same value of DM results in a larger increase in the mortality risk in males compared to females). Increasing patterns of Q(t) in both sexes show that the same value of deviations from the norm (i.e., the same value of DM) causes a larger increase in mortality rate (compared to the baseline rate at respective age) at old ages than it does at younger ages. This is true for both sexes but a faster increase of Q(t) with age in males implies that this “additional price” which an organism has to pay in terms of increasing mortality rate for deviations of physiological parameters from the “norm” increases faster with age in males than in females. This, along with a higher baseline mortality rate, results in the higher mortality rate for males that is observed in human populations.
Figure 1. Estimates of different components of stochastic process model applied to different variants of DM. (A) baseline mortality rate μ0(t); (B) multiplier in quadratic hazard term Q(t); (C) mean allostatic trajectory f1(t); (D) absolute value of feedback coefficient a(t). Different lines correspond to estimates of respective components in specific DM variants for females (solid lines) and males (dashed lines).
Figure 2. Box plots of parameter estimates in different variants of DM applied to data on females (F) and males (M). Note: points are shown as outliers (“+”) if they are larger than q3 + w(q3 − q1) or smaller than q1 − w(q3 − q1), where q1 and q3 are the 25th and 75th percentiles, respectively, and w = 1.5.
Mean allostatic trajectories f1(t) show the opposite pattern – they are higher and increase faster in females. This indicates that the physiological variables summarized by the respective DM tend to deviate farther from the “norm” in females and with age this gap widens faster in females. As such deviations generally correspond to abnormal values of physiological variables that are indicators of diseases/conditions (e.g., diabetes or hypertension), this can contribute to the observed higher prevalence of aging-related diseases and conditions in females. These results on sex differences agree with previous studies showing higher mortality in males but greater propensity for clinical frailty in females (38–40). However, the effect of this on mortality is attenuated by the observed lower values of the multiplier in the quadratic hazard term Q(t) in females compared to males. Nevertheless, for some compositions of DM, we can see no major differences in mean allostatic trajectories for females and males. This observation implies that for some specific indicators of health the male–female mortality–morbidity paradox (41, 42) may not hold and may even be reversed (43).
The results on the age dynamics of the feedback coefficient a(t) are mixed, although females have a tendency to have worse adaptive capacity (both initial adaptive capacity at younger ages and adaptive capacity at old ages). But generally, it confirms our earlier observations that there is no universal behavior of the decline in adaptive capacity for different physiological variables in females and males (4, 26).
Figure 2 displays box plots of parameter estimates in different variants of DM applied to female and male data. This figure illustrates the sex differences in the components of the model described above from the perspective of model’s parameters. It also shows the estimates of parameters σ0 and σ1 indicating that females and males also differ in terms of variability of DM (both at baseline and dynamic variability over age), and, respectively, in terms of variability of the underlying physiological variables, with females having higher variability than males.
Advantages, Limitations, and Further Perspectives
The advantages of the presented mathematical model of age trajectories of physiological dysregulation, aging, and mortality combine those of the underlying model (SPM) and the measure of physiological dysregulation (DM), i.e., a robust measure of physiological dysregulation with a biologically explicit and easily interpretable model of age-related changes in physiology in relation to mortality (4, 5, 8, 10, 26, 27). An additional advantage is that the use of this measure of physiological dysregulation (DM) allows us to perform analysis of a one-dimensional SPM instead of using its multidimensional version.
Physiological variables can be analyzed within the SPM approach by using them explicitly in the hazard or in the summary measures such as the distance DM. The first approach is more flexible in the sense that we have the capability to specify parameters of the model for each variable separately and, thus, to investigate their effects on the time-to-event outcome in detail, as well as to make inference on their dynamic properties from respective parameters in Eq. 3. However, this flexibility has its price as we deal with a multidimensional model and the computational workload in the likelihood estimation procedure is essential. Therefore, if one is interested in the combined effect of the variables on the outcome of interest (and also is not investigating the dynamic properties of each individual variable) then the approach implementing the summary measures such as DM becomes beneficial. It substantially reduces the computational burden as it works in a one-dimensional setting while it still has the same components, which can be interpreted in the context of aging-related processes. In the case of a large number of physiological variables this may result in a considerable computational advantage because the calculations of DM are much faster than the maximization of the likelihood function of the multidimensional SPM for the original set of physiological variables.
From a theoretical standpoint, use of a global measure of physiological dysregulation is likely to be more robust and less noisy than individual biomarkers, which may vary for many reasons unrelated to long-term aging processes. The interpretation of DM as a robust measure is enhanced by the fact that our results here are highly concordant with previous studies on the Women’s Health and Aging Study (WHAS), the Baltimore Longitudinal Study on Aging (BLSA), and InCHIANTI (5, 20, 44), despite major differences in (a) statistical approach; (b) biomarker choice (only three of the biomarkers here – CH, glucose, and hematocrit – were among the 40+ used in those studies); (c) participant characteristics (e.g., WHAS contained exclusively older women from the Baltimore area; InCHIANTI was in Italy); and (d) follow-up time (~5–10 years for WHAS, BLSA, and InCHIANTI versus ~55 years here). This replication of results confirms the idea that DM represents an underlying system-level property of physiological dysregulation that can be detected and interpreted relatively robustly under different contexts and with different markers.
The proposed model provides an approach to measure the impact of systemic dysregulation in an organism on mortality risk, which is an alternative to that using a cumulative measure of health deterioration such as an index of cumulative deficits or deficits index (DI) (45, 46), also an efficient approach to investigate aging-related processes of health deterioration. The DI, which takes into account the cumulative contribution of different variables (including those with possibly minor effects of individual variables on the risk) on mortality risk, can also be implemented in the SPM as in our earlier work (22). Remarkably, both approaches revealed similar sex-specific differences in dynamics of model’s components, see Figure 1 here and Figure 1 in Yashin et al. (22): baseline mortality is higher in males; U-shape (J-shape) of mortality as a function of DI (DM) is narrower and it narrows faster with age in males; and mean allostatic trajectories are higher and increase faster in females (with respective trajectories estimated at zero for males in the case of DI). Such similarities are observed despite, again, major differences in (a) the measure used in the model (DM and DI); (b) variables used in construction of respective measures (biomarkers here versus questions from the survey questionnaires in the DI studies); (c) participant characteristics, sample size, and follow-up time (the model with DI was applied to the National Long Term Care Survey (NLTCS) data, a nationally representative survey of more than 49,000 Medicare enrollees); and (d) specifications of the SPM [gamma-Gompertz (logistic) baseline hazards, non-symmetric U-shapes, non-zero “norms,” quadratic functions for mean allostatic trajectories, and constant feedback coefficients in the SPM applied to the NLTCS data]. This similarity in the results obtained using these two different approaches can thus indicate manifestation of different aspects of the same general process of aging-related deterioration in an organism that cause sexual dimorphism in the dynamics of aging-related processes.
There are two limitations in the implementation of the model described in this paper. First, we note that, as a single imputation method was used, the uncertainty and SEs could in general be underestimated. Another limitation is that it defines the “normal” state of physiological variables from which deviations are evaluated by DM using the mean values at the baseline exam in individuals 40 years and younger. That is, it is assumed that these values minimize the risk of death (as a function of these physiological variables) at all ages. Note that, if we insert the formula of DM (Eq. 1) in Eq. 4 then, essentially, the vector of means in DM represents the “physiological norm” or “optimum” for respective physiological variables in the corresponding multidimensional SPM for the original variables included in the definition of DM. Thus, we actually assume here that this “physiological optimum” is age independent. While this assumption can be true for some variables, for other variables it is likely to be incorrect. Many hormone-related variables experience changes in accordance with the ontogenetic program, and these changes are likely to modify the optimal value for these variables [e.g., this could be manifested by the menopause; see also discussion on age-dependent physiological norms in Yashin et al. (3)]. Thus, one possibility to enhance the results presented here is to assume age-dependent “optima” in the definition of DM. However, this can generally be a challenging problem because one needs to decide which values to use to define the “optima” for different ages before performing analyses of SPM rather than estimate these “physiological optima” from the SPM itself as in our earlier applications (3, 4, 21–24). A possible remedy could be to use a “surrogate” definition of the “physiological optimum” suggested in Yashin et al. (26) taking the average age trajectories of physiological variables for long-lived individuals (say, 90+) on the premise that the long-lived individuals are those who, on the average, managed to keep the age trajectories of physiological variables close to the optimal ones that minimize the risk of death at respective ages. However, there are still computational challenges in this approach as one will need a large enough number of long-lived individuals to reliably estimate the means and the variance–covariance matrices in respective age groups as needed for DM. Alternatively, we can relax the assumption on zero “optimum” f0(t) (and this is also a testable hypothesis in this approach). Although in the current implementation of the model a non-zero “optimum” for DM could mean different things depending on the patterns of changes in the underlying physiological variables, this approach could provide insights on whether the values of biomarkers that are “optimal” for younger individuals are still “optimal” at older ages.
Another direction for applications of the approach presented in this paper is to investigate genetic effects of different candidate genes or single-nucleotide polymorphisms (SNPs) on the measures of physiological dysregulation and on respective aging-related characteristics such as mean allostatic trajectory, resistance to stresses or adaptive capacity which can be evaluated in SPM. For example, one can evaluate carriers of which alleles or genotypes have higher levels of physiological dysregulation, a faster increase in systemic dysregulation with age, a faster decline in resistance to stresses or adaptive capacity calculating respective estimates of SPM components for carriers of different alleles or genotypes. This can be done using either a stratified analysis by genotype or allele using the original SPM (8) as we did here stratifying by sex or implementing DM in the genetic SPM (9), which has an advantage of using additional information on non-genotyped individuals to increase the power. This can also be done either for individual genes/SNPs or for cumulative “genetic doses,” that is the variable calculating the number of some alleles in an individual’s genome that were pre-selected using some criteria or analyses [e.g., “longevity alleles” as in Yashin et al. (27)].
The desire to have useful biological interpretation of the results of analyses stimulates the need for further development of concepts and models describing the aging-related changes developing in the human body. Although the “dysregulation” variable (DM) is useful because it allows us to quantify aging-related changes showing how the particular individual differs from the “ideal” standard (e.g., observed in healthy young adults), its applications to studying mechanisms of aging-related changes requires its further development. From the point of view of this paper, the dysregulation is the result of everything that causes the constructed measure (DM) to deviate from zero where the zero value characterizes the healthy young adults. At the same time, such measure could be too general if one would like to understand causes and mechanisms of aging-related changes, for example, to produce useful recommendations for improving person’s health. One approach might be to examine DM by specific physiological system (47). Alternatively, in Ukraintseva and Yashin (48) and Arbeev et al. (49), the importance of distinguishing among three components of aging-related changes was emphasized. These include basal (senescent), ontogenetic, and exposure related components. The approach described above has a potential to include these components into the model and use them to address questions about various mechanisms driving dynamics of aging-related changes during the life course. For example, allostatic load can be linked with the exposure-related component. Changes in physiological norm can be linked with the ontogenetic component and changes in stress resistance and adaptive capacity with the basal (senescent) component. All three components contribute to dysregulation in physiological state, although their interaction may sometimes have side effects that might slow down or even reduce DM. For example, the reduction of adaptive capacity may slow down accumulation of DM due to weakening negative feedback, especially when f1(t) continues to increase. This process could partly be due to compensatory adaptation. The effects, however, may not last long because in case of weak feedback regulation the “noise” component in Eq. 3 will make stronger effects on Y(t). Note that different components of aging-related changes can be evaluated from available longitudinal data (24, 27). More detailed analyses of dysregulation components will be performed elsewhere.
Implementation of the measure of physiological dysregulation (Mahalanobis distance) in the framework of the SPM of aging allows one to investigate physiological dysregulation in relation to different hidden mechanisms of aging-related changes, and, ultimately, to yield estimates of how such dynamic relationships can produce an increase in the risk of death with age and how it may be related to the observed sex-specific differences in mortality risks. Results of application of the method to individuals from the Framingham original cohort indicate that physiological dysregulation increases with age; that increased dysregulation results in increased mortality, and increasingly so with age; and that, in most but not all cases, there is a decreasing ability to return to baseline physiological state with age. We also show substantial sex differences in these processes, with women becoming dysregulated more quickly but with men showing a much greater sensitivity to dysregulation in terms of mortality risk.
Designed the study: KA and AC; performed statistical analyses and data preparation: KA and LA; contributed software for computation of statistical distance: AC and EM; contributed to methodological work with stochastic process model implementing statistical distance: KA, ES, IA, and AY; wrote the paper: KA, AC, and EM; contributed to the final version and interpretation of results: ES, AK, SU, KC, and AY.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The Framingham Heart Study (FHS) is conducted and supported by the National Heart, Lung, and Blood Institute (NHLBI) in collaboration with Boston University (Contract No. N01-HC-25195). This manuscript was not prepared in collaboration with investigators of the Framingham Heart Study and does not necessarily reflect the opinions or views of the Framingham Heart Study, Boston University, or NHLBI. Research reported in this publication was partly supported by the National Institute on Aging of the National Institutes of Health (NIA/NIH) under Award Numbers R01AG030612, R21AG045245, R01AG046860, P01AG043352, and P30AG034424 (KA, LA, ES, AK, IA, SU, and AY). The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. This work was also partly supported by the NIA/NIH grant 5U01AG023712-10 “Long Life Family Study (LLFS): Duke University/University of Southern Denmark Field Center” (KA, LA, ES, AK, SU, KC, and AY). This paper was not prepared in collaboration with investigators from other LLFS Field Centers, and thus does not necessarily reflect the opinions or views of the entire LLFS group. AC is a member of the FRQ-S-supported Centre de recherche sur le vieillissement and Centre de recherche du CHUS and is a funded Research Scholar of the FRQ-S. This research was additionally supported by CIHR grant #s 110789, 120305, 119485.
The Supplementary Material for this article can be found online at https://www.frontiersin.org/article/10.3389/fpubh.2016.00003
1. Fried LP, Xue Q-L, Cappola AR, Ferrucci L, Chaves P, Varadhan R, et al. Nonlinear multisystem physiological dysregulation associated with frailty in older women: implications for etiology and treatment. J Gerontol A Biol Sci Med Sci (2009) 64:1049–57. doi: 10.1093/gerona/glp076
3. Yashin AI, Ukraintseva SV, Arbeev KG, Akushevich I, Arbeeva LS, Kulminski AM. Maintaining physiological state for exceptional survival: what is the normal level of blood glucose and does it change with age? Mech Ageing Dev (2009) 130:611–8. doi:10.1016/j.mad.2009.07.004
4. Arbeev KG, Ukraintseva SV, Akushevich I, Kulminski AM, Arbeeva LS, Akushevich L, et al. Age trajectories of physiological indices in relation to healthy life course. Mech Ageing Dev (2011) 132:93–102. doi:10.1016/j.mad.2011.01.001
5. Cohen AA, Milot E, Yong J, Seplaki CL, Fueloep T, Bandeen-Roche K, et al. A novel statistical approach shows evidence for multi-system physiological dysregulation during aging. Mech Ageing Dev (2013) 134:110–7. doi:10.1016/j.mad.2013.01.004
7. Cohen AA, Milot E, Li Q, Bergeron P, Poirier R, Dusseault-Belanger F, et al. Detection of a novel, integrative aging process suggests complex physiological integration. PLoS One (2015) 10(3):e0116489. doi:10.1371/journal.pone.0116489
8. Yashin AI, Arbeev KG, Akushevich I, Kulminski A, Akushevich L, Ukraintseva SV. Stochastic model for analysis of longitudinal data on aging and mortality. Math Biosci (2007) 208:538–51. doi:10.1016/j.mbs.2006.11.006
9. Arbeev KG, Akushevich I, Kulminski AM, Arbeeva LS, Akushevich L, Ukraintseva SV, et al. Genetic model for longitudinal studies of aging, health, and longevity and its potential application to incomplete data. J Theor Biol (2009) 258:103–11. doi:10.1016/j.jtbi.2009.01.023
10. Yashin AI, Arbeev KG, Akushevich I, Kulminski A, Ukraintseva SV, Stallard E, et al. The quadratic hazard model for analyzing longitudinal data on aging, health, and the life span. Phys Life Rev (2012) 9:177–88. doi:10.1016/j.plrev.2012.05.002
18. Singer BH, Ryff CD, Seeman T. Operationalizing allostatic load. In: Schulkin J, editor. Allostasis, Homeostasis, and the Costs of Physiological Adaptation. Cambridge: Cambridge University Press (2004). p. 113–49.
19. Milot E, Cohen AA, Vezina F, Buehler DM, Matson KD, Piersma T. A novel integrative method for measuring body condition in ecological studies based on physiological dysregulation. Methods Ecol Evol (2014) 5:146–55. doi:10.1111/2041-210X.12145
20. Milot E, Morissette-Thomas V, Li Q, Fried LP, Ferrucci L, Cohen AA. Trajectories of physiological dysregulation predicts mortality and health outcomes in a consistent manner across three populations. Mech Ageing Dev (2014) 141:56–63. doi:10.1016/j.mad.2014.10.001
21. Yashin AI, Arbeev KG, Kulminski A, Akushevich I, Akushevich L, Ukraintseva SV. Cumulative index of elderly disorders and its dynamic contribution to mortality and longevity. Rejuvenation Res (2007) 10:75–86. doi:10.1089/rej.2006.0500
22. Yashin AI, Arbeev KG, Kulminski A, Akushevich I, Akushevich L, Ukraintseva SV. Health decline, aging and mortality: how are they related? Biogerontology (2007) 8:291–302. doi:10.1007/s10522-006-9073-3
23. Yashin AI, Arbeev KG, Kulminski A, Akushevich I, Akushevich L, Ukraintseva SV. What age trajectories of cumulative deficits and medical costs tell us about individual aging and mortality risk: findings from the NLTCS-medicare data. Mech Ageing Dev (2008) 129:191–200. doi:10.1016/j.mad.2007.12.005
24. Yashin AI, Arbeev KG, Akushevich I, Ukraintseva SV, Kulminski A, Arbeeva LS, et al. Exceptional survivors have lower age trajectories of blood glucose: lessons from longitudinal data. Biogerontology (2010) 11:257–65. doi:10.1007/s10522-009-9243-1
25. Arbeev KG, Ukraintseva SV, Kulminski AM, Akushevich I, Arbeeva LS, Culminskaya IV, et al. Effect of the APOE polymorphism and age trajectories of physiological variables on mortality: application of genetic stochastic process model of aging. Scientifica (Cairo) (2012) 2012:568628. doi:10.6064/2012/568628
26. Yashin AI, Arbeev KG, Ukraintseva SV, Akushevich I, Kulminski A. Patterns of aging related changes on the way to 100: an approach to studying aging, mortality, and longevity from longitudinal data. N Am Actuar J (2012) 16:403–33. doi:10.1080/10920277.2012.10597640
27. Yashin AI, Arbeev KG, Wu D, Arbeeva LS, Kulminski A, Akushevich I, et al. How lifespan associated genes modulate aging changes: lessons from analysis of longitudinal data. Front Genetics (2013) 4:3. doi:10.3389/fgene.2013.00003
31. Yashin AI, Manton KG, Stallard E. The propagation of uncertainty in human mortality processes operating in stochastic environments. Theor Popul Biol (1989) 35:119–41. doi:10.1016/0040-5809(89)90013-0
32. Yashin AI. Dynamics in survival analysis: conditional Gaussian property vs. Cameron-Martin formula. In: Krylov NV, Lipster RS, Novikov AA, editors. Statistics and Control of Stochastic Processes. New York: Springer (1985). p. 446–75.
34. Cohen AA, Li Q, Milot E, Leroux M, Faucher S, Morissette-Thomas V, et al. Statistical distance as a measure of physiological dysregulation is largely robust to variation in its biomarker composition. PLoS One (2015) 10(4):e0122541. doi:10.1371/journal.pone.0122541
35. Yashin AI, Arbeev KG, Akushevich I, Arbeeva L, Kravchenko J, Il’yasova D, et al. Dynamic determinants of longevity and exceptional health. Curr Gerontol Geriatr Res (2010) 2010:381637. doi:10.1155/2010/381637
36. Manton KG, Stallard E, Woodbury MA, Dowd JE. Time-varying covariates in models of human mortality and aging: multidimensional generalizations of the Gompertz. J Gerontol (1994) 49:B169–90. doi:10.1093/geronj/49.4.B169
38. Mitnitski A, Song X, Skoog I, Broe GA, Cox JL, Grunfeld E, et al. Relative fitness and frailty of elderly men and women in developed countries and their relationship with mortality. J Am Geriatr Soc (2005) 53:2184–9. doi:10.1111/j.1532-5415.2005.00506.x
39. Puts MT, Lips P, Deeg DJ. Sex differences in the risk of frailty for mortality independent of disability and chronic diseases. J Am Geriatr Soc (2005) 53:40–7. doi:10.1111/j.1532-5415.2005.53008.x
43. Kulminski AM, Culminskaya IV, Ukraintseva SV, Arbeev KG, Land KC, Yashin AI. Sex-specific health deterioration and mortality: the morbidity-mortality paradox over age and time. Exp Gerontol (2008) 43:1052–7. doi:10.1016/j.exger.2008.09.007
44. Cohen AA, Milot E, Li Q, Legault V, Fried LP, Ferrucci L. Cross-population validation of statistical distance as a measure of physiological dysregulation during aging. Exp Gerontol (2014) 57:203–10. doi:10.1016/j.exger.2014.04.016
46. Kulminski A, Yashin A, Arbeev K, Akushevich I, Ukraintseva S, Land K, et al. Cumulative index of health disorders as an indicator of aging-associated processes in the elderly: results from analyses of the national long term care survey. Mech Ageing Dev (2007) 128:250–8. doi:10.1016/j.mad.2006.12.004
47. Li Q, Wang S, Milot E, Bergeron P, Ferrucci L, Fried LP, et al. Homeostatic dysregulation proceeds in parallel in multiple physiological systems. Aging Cell (2015) 14:1103–12. doi:10.1111/acel.12402
Keywords: physiological dysregulation, stochastic process model, Mahalanobis distance, longitudinal data, mortality, sex differences
Citation: Arbeev KG, Cohen AA, Arbeeva LS, Milot E, Stallard E, Kulminski AM, Akushevich I, Ukraintseva SV, Christensen K and Yashin AI (2016) Optimal Versus Realized Trajectories of Physiological Dysregulation in Aging and Their Relation to Sex-Specific Mortality Risk. Front. Public Health 4:3. doi: 10.3389/fpubh.2016.00003
Received: 04 December 2015; Accepted: 11 January 2016;
Published: 25 January 2016
Edited by:Onyebuchi A. Arah, University of California Los Angeles (UCLA), USA
Reviewed by:Arnold Mitnitski, Dalhousie University, Canada
Catherine M. Crespi, University of California Los Angeles (UCLA), USA
Copyright: © 2016 Arbeev, Cohen, Arbeeva, Milot, Stallard, Kulminski, Akushevich, Ukraintseva, Christensen and Yashin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Konstantin G. Arbeev, firstname.lastname@example.org