Optimal Versus Realized Trajectories of Physiological Dysregulation in Aging and Their Relation to Sex-Specific Mortality Risk

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.


inTrODUcTiOn
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)(6)(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)(9)(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)(12)(13)(14)(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 twodimensional cases (3,4,(21)(22)(23)(24)(25)(26)(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 onedimensional 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.

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 X 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 agedependence 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 variancecovariance 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  (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: 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: ln ( ) ln ( ) σ . 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 agingrelated 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)  (i.e., b f 1 0 = ). 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. " 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: f t a f 1 1 ( ) = 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 b f 1 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   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. Mean allostatic trajectories f1(t) show the opposite patternthey 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 agingrelated 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)(39)(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 malefemale 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 X 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)(22)(23)(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 longlived 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.
cOnclUsiOn 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.

aUThOr cOnTriBUTiOns
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,