Network dynamical stability analysis reveals key “mallostatic” natural variables that erode homeostasis and drive age-related decline of health

Using longitudinal study data, we dynamically model how aging affects homeostasis in both mice and humans. We operationalize homeostasis as a multivariate mean-reverting stochastic process. We hypothesize that biomarkers have stable equilibrium values, but that deviations from equilibrium of each biomarker affects other biomarkers through an interaction network—this precludes univariate analysis. We therefore looked for age-related changes to homeostasis using dynamic network stability analysis, which transforms observed biomarker data into independent “natural” variables and determines their associated recovery rates. Most natural variables remained near equilibrium and were essentially constant in time. A small number of natural variables were unable to equilibrate due to a gradual drift with age in their homeostatic equilibrium, i.e. allostasis. This drift caused them to accumulate over the lifespan course and makes them natural aging variables. Their rate of accumulation was correlated with risk of adverse outcomes: death or dementia onset. We call this tendency for aging organisms to drift towards an equilibrium position of ever-worsening health “mallostasis”. We demonstrate that the effects of mallostasis on observed biomarkers are spread out through the interaction network. This could provide a redundancy mechanism to preserve functioning until multi-system dysfunction emerges at advanced ages.


Introduction
Homeostasis is the self-regulating process that maintains internal stability 1 .Yet as individuals age, it is characteristic for biomarkers to drift away from healthy levels; something about homeostasis is therefore "lost" during the aging process 2 .For example, loss of protein homeostasis is believed to cause the hallmark accumulation of unfolded, misfolded and aggregate proteins with age 3 .Accumulation is observed at multiple biological scales, including oxidative damage 4 , epigenetic age 5 , senescent cells 6 , and regulatory T-cells 7 at the cellular scale, and extending up to the whole organism scale where clinical deficits 8 , including chronic diseases 9 , accumulate with age.Sehl and Yates performed univariate analysis of 445 health biomarkers and found that almost all of them accumulate negatively with age -typically showing linear decline 10 .Such accumulation of biomarker values in a particular direction appears to be a generic feature of aging.When biomarkers reach abnormal values, they are associated with dysfunction and poor health, independently of age 11,12 .A general mechanism of how accumulation and poor health emerge from homeostasis has, nevertheless, been missing.
Prior work suggests that accumulation may be a consequence of a drifting equilibrium position.Allostasis, literally "homeostasis through change" 13 , describes a version of homeostasis in which the equilibrium position is mutable, adapting as necessary to environmental demands 14 .Over time, "wear-and-tear" of this adaptive stress-response causes a subclinical accumulation of dysfunction known as "allostatic load" 14 .We hypothesize that these allostatic changes may be asymmetric, causing a coherent, population-level drift in equilibrium biomarker values with age, and ultimately leading to accumulating biomarker values in particular directions.
Directly estimating an individual's allostatic load remains an open challenge 14 , owing to the confounding effects of the underlying interaction networks 12 .Instead, most algorithms infer allostatic load by outlier detection 12,14 or other symmetric indicators, agnostic to any preferred biomarker accumulation direction 15,16 .These approaches have not been reconciled with generic, age-associated biomarker accumulation, which proceeds in preferred directions 10 .It therefore remains unclear how allostatic load leads to worsening health.Other theories posit that outlying biomarker values indicate damage, which promotes further damage e.g. as quantified by the number of health deficits ("frailty index", FI) 11,17,18 .Needed is direct evidence of allostasis and how it is associated with worsening health.
Instability is another mechanism for accumulation.While linear accumulation is the norm 10 , some biomarkers accumulate exponentially with age e.g.senescent cells 6 and the FI 8 .Exponential growth indicates an instability 19 .Nevertheless, a weak arXiv:2306.12448v2[q-bio.OT] 23 Oct 2023 instability can appear linear until advanced ages.As a result, it remains unclear whether age-related accumulation proceeds due to a shifting equilibrium, a weak instability, or some hybrid of the two.
Operationalizing and quantifying homeostatic changes is challenging 14 because homeostasis is a property of the whole system, not individual constituent parts 1,12 .In the language of complexity science 20 , homeostasis is an emergent property of a network of interacting variables.Each variable measures a part of the system, but changes to one part can be balanced by other parts.For example, heart rate declines with age but can be compensated for by increased stroke volume 10 in order to maintain arterial blood pressure 1 .The essential aspects of homeostasis are: (i) a multivariate interacting dynamical system, (ii) an equilibrium state, which may vary with age (allostasis), (iii) the system spontaneously returns to the equilibrium state (dynamical stability), and (iv) stresses (and interventions) provide random shocks to the system.Altogether, homeostasis can be operationalized as a multivariate, mean-reverting stochastic process 16 .
Dynamical stability analysis uses eigen-decomposition to probe the stability of arbitrary systems 21,22 .The system is first linearized around an equilibrium position 21 .Orthogonal eigenvectors are then identified that decouple the interactions between variables.Eigenvectors are composite health measures that serve as natural variables since they do not interact or compensate for each other, and so can be analyzed individually.Each such natural variable has an associated eigenvalue that determines its stability via a characteristic recovery rate or timescale (−eigenvalue = rate = timescale −1 ).A system is stable if and only if all recovery rates are positive 21 .Conversely, dynamical instability arises only if at least one recovery rate is negative.
We confront homeostasis with minimal assumptions.We seek generic changes to biomarker equilibrium and stability within aging organisms.We investigate multiple longitudinal datasets with multiple organisms (mice and humans) and multiple outcomes (dementia and death).In contrast to earlier work by Sehl and Yates 10 , our model is multivariate and generic so that we can model homeostasis without constraining its dynamical behaviour.We find that allostatic drift is consistent with the observed data.Importantly, we find that a small set of natural variables drive mortality and can be used to characterize an individual's health state.We do not observe any dynamical instabilities.

Model
To analyse stability for deterministic 21 or stochastic 22 dynamics, we use a linear approximation near a stable point, ⃗ µ in ≡ ⃗ µ 0 + Λ Λ Λ⃗ x in +⃗ µ age t in (1)   where i indexes the individual, n indexes the timepoint, t n is the age, ⃗ y is a vector of observed biomarkers, and ⃗ x is a vector of covariates that includes sex.W W W , Σ Σ Σ, and Λ Λ Λ are constant matrices, independent of i and n.If we take the average over individuals (indicated by angled brackets) then we can obtain rates of average change We see that changes to the mean of a particular biomarker, y j , are due either to recovery of y j towards the equilibrium position, µ j , or because of interactions with a compensating variable, y k̸ = j through off-diagonal elements of W W W .This provides both a mechanism for biological redundancy -if the organism can actively influence some of the y k then it can use them to steer others -and a mechanism for mutual dysfunction -since W W W couples dysregulation of y k̸ = j to that of y j .In Supplemental Section S8 we show that equation (1) approximates general nonlinear dynamics 21 .We also specifically show that it approximates the stochastic process model 16 , a framework for aging biomarker dynamics.Note that the model permits unequally-spaced sampling of individuals through ∆t in+1 , which is the time between measurements of individual i at time t n and t n+1 .
The stability of the model depends on the eigenvalues of W W W .We can decouple variable means with the eigenvector transformation matrix P P P. We obtain z i jn+1 = z i jn + λ j ∆t in+1 (z i jn − μi jn ) + εi j , where⃗ z n ≡ P P P −1 ⃗ y n , λ j ≡ P −1 j• W W W P • j , ⃗ µ n ≡ P P P −1 µ n and ⃗ ε ≡ P P P −1 ⃗ ε.We refer to⃗ z as natural variables.The natural variables build correlations only through the noise term, ε -in addition to any correlated initial conditions.The system is mutually-diagonal if ε is uncorrelated.While our dynamics are discrete, it is also helpful to consider continuous dynamics corresponding to the limit ∆t → 0; see Figure 1 and Box 1 (also Supplemental Section S8 for more details).
The parameters W W W , ⃗ µ 0 , ⃗ µ age and Λ Λ Λ are estimated from the data (Σ Σ Σ can also be).The stochastic term,⃗ ε, is assumed to be normally distributed and independent across timepoints.See Supplemental Section S6 for details.(For the remainder of the paper, we simplify notation by dropping the tilde and suppressing the individual i and timepoint n indices.)Optimal parameter values are selected by maximizing the likelihood.For uncorrelated noise this reduces to weighted linear regression.
If a mutually-diagonal system reaches steady-state -having run long enough to forget initial conditions -then the natural variables, ⃗ z, are the principal components, ranked by stability (Supplemental equation (S49)).We used principal component analysis (PCA) as a preprocessing step.If P P P is orthogonal (which it is from PCA) then Parseval's theorem states that ⟨∑ j y 2 jn ⟩ = ⟨∑ j z 2 jn ⟩ = ∑ j (Var(z j ) + ⟨z j ⟩ 2 ): this means that a single z k with large mean and variance can dominate that of the ⃗ y.Consider a 1-dimensional space, z.If we take the limit ∆t → 0 then equation ( 3) is a modified Ornstein-Uhlenbeck process.The mean and variance are solutions to ordinary differential equations.The mean is described by where µ age t is the time-dependent part of µ n and µ 0 is the remaining part.The general solution of equation ( 4) is where ⟨z 0 ⟩ is the initially observed mean at t = 0.The exponential factors dampen or aggregate the mean depending on the sign of λ .If λ < 0 the system is stable and once |λ |t ≫ 1 a dynamical steady-state (ss) is reached, The steady-state is equivalent to the system forgetting its initial conditions.This steady-state behaviour can explain the drift observed by Sehl and Yates 10 (Supplemental Section S8.6).In the steady-state, the mean drifts at a constant rate, but there is a constant lag of ⟨z − µ⟩ ss = µ age /λ ; Figure 1 illustrates.Only when µ age = 0 (no drift) is µ(t) the steady-state position.Outside of steady-state, the mean is displaced by for reference time t 0 .The first term encodes the system's initial conditions, whereas the last term encodes long-time drifting behaviour.Systems near steady-state exhibit Memory ≪ Drift.If λ = 0 the system is marginally stable and preserves its initial conditions.If λ > 0 the system is unstable and the initial conditions grow exponentially over time.In either case the steady-state is never reached.In contrast to the mean, for λ < 0 the variance eventually equilibrates, reaching a constant value.The variance is described by where σ 2 is the noise strength.The general solution is given by Approaching instability, with λ → 0, the system accumulates noise.

Results
We analysed four datasets originating from three studies: two human and two mouse.Our analysis focused on the key properties of homeostasis: stability and equilibrium position.We used model selection to compare our model to the null hypothesis and to pick an optimal model form (Supplemental Section S5).
We observed that both an interaction network, W W W , and an equilibrium term, ⃗ µ, were needed to optimally predict future biomarker values.We saw no evidence of nonlinear terms in the dynamics.We found that fitting equation (3) using principal components (PCs) yielded equivalent performance to the model with full flexibility, equation ( 1), but was already in diagonal form.Hence for each dataset we analysed a set of decoupled, one-dimensional equations in z j (with j sorted by stability, so that j = 1 is the least stable in each study).
For covariates, we generally found non-significant improvements in prediction -though we kept them to improve interpretability (to reduce confounding).The exception was the age covariate, µ age , which significantly improved the fit of the SLAM Het3 mice (SLAM C57/BL6 were almost significant).The presence of µ age indicates allostasis in the form of a time-dependent homeostasis.
The interaction networks between variables can be represented by the respective weight matrices, e.g. Figure 2A.For ELSA we see expected relationships, e.g.total/HDL/LDL cholesterol, and non-dominant/dominant grip strength.ELSA also shows a block of like-variables including the FI-ADL, FI-IADL, self-reported health (srh), and gait speed, which could relate to frailty 23 .See Supplemental Figure S9 for networks from the other datasets.Recovery rates in human-equivalent (h.e.) years i.e. negative eigenvalues (−λ ).The smallest recovery rates determine system stability 21 .A recovery rate of 0.025 implies 1 − e −1 = 63% recovery after −λ −1 = 40 years (95% recovery after 120 years).
The survival data all have similar minimum rates near 0.025, whereas the dementia data was faster (Paquid).The dotted lines are network diagonals (−W j j ); the solid lines are rates (−λ j ).
The interactions between observed biomarkers prevents us from assessing stability directly.However, we can eigendecompose the networks to yield an equivalent non-interacting network of natural variables.Each natural variable has a characteristic recovery rate, Figure 2B.All natural variables were stable, with λ < 0. Faster recovery rates indicate higher stability (resilience) 22 .It takes 3 timescales for the system mean to recover 95% of the way to equilibrium.For each mortality dataset the recovery timescale of its slowest natural variable was comparable to the organism's lifespan, ≈ 40 human-equivalent years; only the mental acuity dataset (Paquid) was faster (≈ 20 years).In all datasets the rates for the natural variables extended to higher and lower values than the diagonal elements of the observed biomarkers (compare solid to dotted lines) -indicating that network interactions play an important role in recovery dynamics.Position relative to equilibrium vs recover rate.Most natural variables were homeostatic (near equilibrium at 0).Some (labeled) variables were observed to be far from equilibrium; variables are labelled by rank e.g.01 ≡ z 01 has the fastest recovery (furthest left).B. Characterization of natural variable deviations from equilibrium using equation (8).Observe that ELSA is the only dataset where memory may dominate the system behaviour (ratio ≲ 1 = 10 0 ), indicating that the followup period may have been too short to reach a steady-state.In both figures only mouse (SLAM) data points over age 80 weeks were used since biomarkers had a u-shaped curve over the lifespan 24 .
We summarize homeostasis in Figure 3A, using the population means.If variables are in homeostasis then the mean should be close to µ n , where the scale is determined by the native dispersion.Each dataset had most natural variables near 0 with a few outliers, such as z 1 for all datasets.(In contrast, the majority of observed biomarkers had large deviations from equilibriumsee Supplemental Figure S10).We characterized the natural variable dynamics using equation (8) in Figure 3B.Excluding ELSA, most data points were in a steady-state, as indicated by their small memory term (relative to drift).The steady-state mean includes a drift caused by µ age .Across variables, the deviations from equilibrium observed in Figure 3A, ⟨z n − µ n ⟩, were very strongly correlated with µ age , with correlation coefficients: -0.988 (p = 2 • 10 −4 , SLAM BL/6), -0.947 (p = 10 −3 , SLAM Het3), -0.989 (p = 0.01, Paquid), and -0.302 (p = 0.14, ELSA).This is consistent with equation ( 6), and supports our use of an allostatic model with equilibrium drifts given by µ age .The smaller correlations observed with the ELSA dataset are consistent with the strong memory effect seen in Figure 3B -violating the steady-state assumption of equation (6).ELSA may have failed to reach steady-state due to the limited followup period, which was the shortest of all datasets by a factor of 2, or could indicate the confounding effects of medical interventions, which are not relevant for the other datasets.
Most natural variables have small drift and are effectively homeostatic -with only a few strongly drifting allostatic natural variables.The steady-state drift rate of natural variables, µ age , was correlated with the survival risk for each dimension: Figure 4A.The correlations were typically strong: -0.958 (p = 0.002, SLAM BL/6), -0.713 (p = 0.1 SLAM Het3), -0.987 (p = 0.01, Paquid), and -0.534 (p = 0.006 ELSA); overall: -0.742 (p = 3 • 10 −8 ).The correlation was weakest for ELSA, which had not reached steady-state.The Cox proportional hazards coefficients, conditioned on age and sex, showed a similarly strong correlation with µ age , 0.70 (p = 10 −7 , all data) (Supplemental Figure S14).Furthermore, we see that the drift direction, sign(µ age ), is the same as the risk direction (p = 0.0003, Fisher test).Hence, not only does homeostasis drift with age, the direction of the drift is towards ill-health.The primary risk directions were z 1 for ELSA and Paquid and z 2 for SLAM.Interestingly, z 2 of the Het3 mice is nearly identical to z 1 of the C57BL/6 mice in terms of covariates and survival effecthence z 1 of the C57BL/6 is also likely a key risk direction (Supplemental Figures S11 and S18).Regardless, z 1 or z 2 exhibited the strongest survival effect for their each dataset (Figure 4A).These variables also both had small eigenvalues (z 1 is rank 1 and z 2 is rank 2).However, this relationship between survival and eigenvalue magnitude does not appear to generalize, see Figure 4B.
As an illustration of the utility of the correlation between survival and µ age , we consider a simple summary health measure.The Cox proportional hazards model assumes the hazard scales as exp ( ⃗ β T ⃗ z), where the jth coefficient, β j , is the log-hazard ratio per unit increase of z j .As mentioned in the previous paragraph, ⃗ β ∼ ⃗ µ age were correlated; the relative hazard can therefore be approximated by exp (⃗ µ T age ⃗ z).Indeed, we observed that b ≡ ⃗ µ T age ⃗ z is an excellent predictor of survival, see Figure 5A and Supplemental Figure S15.
The natural variables with large |µ age | will eventually experience the largest drift, according to equation (7).z 1 in Figure 5B is an example of such an accumulating variable.The other variables with large |µ age | experienced a similar accumulation (Supplemental Figure S13).For an orthogonal transformation such as P P P −1 , the sum of the variance and squared mean is conserved (Parseval's theorem).Natural variables with large means and variances will therefore disproportionately affect the means and variances of observed biomarkers.The effect is demonstrated for ELSA in Figure 5B.As the dominant natural variable drifts it influences observable biomarkers to drift as well.
Slower recovery rates (eigenvalues) take longer to forget perturbations, causing the associated natural variables to accumulate variance due to noise.Recall that the slowest recovery rates were on the order of a lifespan (Figure 2B).The Pearson correlations between the estimated variance and rate (-eigenvalue) were strong: -0.852 (p = 0.03, SLAM BL/6), -0.802 (p = 0.05 SLAM Het3), -0.998 (p = 0.002, Paquid), and -0.764 (p = 9 • 10 −6 ELSA) (log-log scale; see Supplemental Figure S16).Hence the variances we observe at old ages will be dominated by the variables with the smallest eigenvalues, λ (e.g.z 1 and z 2 ).As we have seen before, these variables are often -but not always -strongly associated with adverse effects, depending on the drift rate µ age .This suggests that most of the age-related changes to health were concentrated in a few z k which drive both biomarker drift (mean) and dispersion (variance).Growing variance along these dimensions may capture individual accumulation of stochastic damage, such as genetic damage or disease.

Discussion
We fit a homeostasis model of equilibrium and stability to four longitudinal aging datasets (two mouse and two human) using generic health biomarkers.Our model is lightweight, can be estimated using standard statistical algorithms, and is sufficient to capture essential information about the aging process.Health biomarkers have an equilibrium position, ⃗ µ.Their corresponding stochastic term, which has covariance Σ Σ Σ, represents random stresses that drive individuals away from equilibrium; as well as residual effects such as individual variability and nonlinearities.An interaction network, W W W , pulls individuals towards equilibrium either through recovery (diagonal terms) or by one variable compensating for another (off-diagonal terms).By eigen-decomposing W W W we transformed the dataset into non-interacting, natural variables -which are linear combinations of  S15 for the other datasets.B. Natural variables can drive changes in observable biomarkers.The z 1 mean is accumulating in the negative direction.This accumulation is mapped into observable variables with ⟨P j1 z 1 ⟩ for indicated timepoints each separated by approximately 4 years.The drift direction is overwhelmingly unhealthy: increased disability measures (srh, eye, hear, FI.ADL and FI.IADLhigh is bad), decreased physical ability scores (gait and grip), increased inflammation (crp), increased glucose, etc.The effect of the drift is concentrated in z 1 but dilute across its covariates, which could make the effect of unhealthy z 1 subclinical in the observed biomarkers.All variables are on standardized scale.Similar effects were observed for the other datasets (Supplemental Figure S13).the input biomarkers.This increases interpretability and simplifies analysis.The stability of the system is described by the recovery rates of the natural variables -which are the corresponding eigenvalues with flipped signs, −λ .
We modelled homeostasis as stability around an equilibrium mean value.Stability can only be assessed using the natural variables because the original biomarkers have interactions between them.Homeostasis was violated by some natural variables (Figure 3A).Although most natural variables had average values near the homeostatic equilibrium -indicative of homeostasis -several were far away.We determined that this latter group were out of homeostasis because they were chasing drifting equilibrium positions from behind.This equilibrium drift with age represents allostasis, i.e. a changeable equilibrium position.
Allostatic variables accumulated over the course of the study period as they chased after the drifting equilibrium, systematically increasing or decreasing.This was facilitated by an age-dependent equilibrium position, governed by µ age , and was typically accompanied by a small eigenvalue.The gap between the population and allostatic equilibrium position is governed by µ age /λ , so a small λ enables a large gap such that the entire population drifts coherently towards the moving equilibrium -causing population-level accumulation.This makes the linear drift term, µ age , the primary culprit for causing biomarkers to drift with age.Presumably µ age arises from either (a) the effects of unknown biomarkers/mechanisms not included in the model (i.e.µ age t ≈ ∑ k̸ = j W jk ⟨y ikn − µ ikn ⟩ in equation ( 2) for a set of unknown y ikn ), or (b) asymmetric stressors, which cannot be captured by our symmetric stochastic term (e.g.there is no such thing as negative damage so health deficits skew positive 8 ).
The transformation to natural variables effectively compressed the drifting (accumulating) mean of many variables into a small number of natural variables.The natural variables can be thought of as the underlying cause of the observed biomarker drift (Figure 5B).In this manner, the widely observed age-related decline in biomarkers 10 are governed by a few natural variables -which are not directly observed.The effect is spread out by the transformation, potentially hiding the observed biomarker decline below diagnostic thresholds.This may be a redundancy mechanism: the network permits the biological system to spread out the age-related decline to keep biomarkers in healthy ranges for longer.The trade-off may be that many biomarkers would reach unhealthy ranges concurrently, leading to multisystem dysfunction.For example, the effects of chronic kidney disease are mild and non-specific until the patient nears kidney failure -at which point multisystem failure is imminent, typically leading to death via cardiovascular disease 25 .This tradeoff assumes that diagnostic thresholds represent critical values beyond which deterioration of a biological system accelerates.Univariate dynamical modelling of senescent cell count in mice 6 and E. coli membrane integrity 26 supports the existence of such critical values, where repair mechanisms saturate and decline accelerates.From this perspective, an individual's robustness would depend on their buffer space available to absorb new insults, which could be quantified by the natural variable scores together with the stressor effect strength which should be proportional to the noise σ .
Consistent with this perspective 6 , the allostatic drift rate, µ age , strongly correlated with the mortality/dementia risk associated with each natural variable (Figure 4A and Supplemental Figure S14).Since µ age is the steady-state drift rate of the mean, the steady-state behaviour is continually worsening health due to the drifting mean.Prior work on operationalizing allostasis has neglected the existence of a preferred risk direction, instead using the absolute distance from allostasis as a mortality factor 12,15,16 , irrespective of whether biomarkers are high or low.In contrast, our results indicate that numerous natural variables do, in fact, have preferred risk directions.Aging researchers should be aware of this symmetry breaking.This means that the adaptive changes due to allostasis at best mitigate declining health and, at worst, lead to a further decline in health.We refer to this phenomenon as "mallostasis": the tendency of an aging biological system towards an ever-worsening equilibrium position.We have used this phenomenon both to identify important survival variables and to generate a novel composite health measure.
Our key quantitative results coincide with three key qualitative predictions made by allostatic load theory: (i) a shifting equilibrium position for biomarkers indicative of adaptive changes (allostasis, Figure 3A), (ii) the shift is associated with adverse outcomes (mallostasis, Figure 4A), and (iii) the shift is subclinical due to compensating mechanisms between biomarkers (transformation, Figure 5B) 14 .This is compelling evidence that allostasis is a generic aging phenomenon, rather than being specific to neuroendocrinology.Our proposed composite health measure is therefore a novel estimator of allostatic load.In contrast to conventional estimators 14 , we were able to estimate allostatic drift directly as µ age .Our results rely on using natural variables, which are canonical coordinates that greatly simplify analysis.
Allostatic load is believed to arise from the long-term costs of short-term protection against stressors 13 , making it an example of antagonistic pleiotropy 27 .Alternatively, long-term costs could reflect imperfect repair.Regardless, long-term costs that accumulate in a given direction would lead to the allostatic drift which we have observed and characterized.Furthermore, we observed slow dynamical rates for the dominant mortality-risk natural variables (Figure 2B).Accordingly, the dynamical timescale of these effects are comparable to the organismal lifetime -consistent with long-term costs.
Interestingly, we did not observe any instabilities or nonlinearities.We had expected that "allostatic overload" -the final state of allostatic load theory 14 -would operationalize as an instability.However, instabilities, and their associated exponential growth, are rare among health biomarkers 10 ; although they are observed in summary measures of health such as the FI (frailty index) 8 .Other unstable, FI-like variables can be extracted from generic biomarkers using nonlinear techniques such as a deep neural network 19 , diagnostic thresholds 11 , or quantile-based preprocessing 18,28 .Since we did not see evidence of instabilities or other nonlinearities in the natural variables, the nonlinear embedding or discretization should be considered as a possible cause for observed FI-like instabilities.It may be that biological systems naturally suppress nonlinear effects of aging -obscuring the effects -or, conversely, that aging is primarily a linear phenomenon that slowly pushes individuals towards nonlinear tolerance thresholds for dysfunction/damage, e.g.saturation of repair 6,26 and/or emergence of chronic disease.A non-trivial issue is that exponential growth often appears linear, for example the FI in mice and younger humans (≲ 85 years old) 29 .Nonlinear effects in biomarker dynamics may require special populations, such as the ill or exceptionally old, to be observed.
The key model variables, z 1 and z 2 , dominate the aging process.These natural variables with smaller λ carried the majority of the variance and become the dominant principal components in the steady-state model (Supplemental Figure S16 and equation (S49), respectively).Applying Parseval's theorem, these variables will control the variance of directly observed biomarkers.Since they also dominate the means via allostatic drift, they will determine the aging phenotype that we observe.Both effects get stronger with age.This means that the empirically observed age-related changes in the mean and variance of biomarkers will be predominantly caused by only a few key natural variables.Hence the nearly-universal linear decline in health biomarkers observed by Sehl and Yates 10 may simply be a few declining natural variables spread across the observed biomarkers (Supplemental Section S8.6).Furthermore, this implies that a single dimensional decline can drive many observed biomarkers, which is the foundational assumption of "biological age" estimators 30,31 .Our results provide much needed support for such low dimensional representations of aging -which should become increasingly accurate with advancing age since the means of the key natural variables grow fastest, and their variances grow largest.
The natural variables, z, should be good choices for targeting and monitoring interventions.They are prospective biomarkers with the convenient property that if you can intervene on one it will not affect the others.In contrast, we know from the network of interactions that intervening on any single biomarker is likely to affect many other biomarkers.In the steady-state, the mallostatic drift rate, controlled by µ age , is a proxy for the hazard and therefore identifies the most important targets of intervention.The coefficients of the transformation, P P P, provides both hints at what mechanisms each z j is capturing as well as a map for which biomarkers will be affected by interventions on z j .For example, z 1 of ELSA shares many features with frailty: strong age dependence, large effect in gait, weakness (grip strength), disability and self-reported health, and large survival effect 23 .z 1 is thus a prospective biomarker of frailty and can be used both to monitor an individual's frailty and to engineer interventions.The strong signals we see in Figure 5B for gait, grip strength and activities of daily living are hints that loss of physical fitness is one mechanism by which frailty proceeds and therefore one mechanism by which we can intervene, consistent with a meta-analysis which has shown that physical activity can reduce frailty in humans 32 .Remarkably, we observed other prospective targets in addition to z 1 of ELSA.Given an organism and set of biomarkers, each z j with substantial drift should be considered a prospective intervention target, with the faster drifting being the most important.
We note a few limitations of our study.We assumed linear, time-invariant interactions through the network, W W W -following previous work that suggested that interactions are linear and time-invariant 33 (as are the principal components 23 ).The networks we extracted were symmetric and hence acausal due to our use of PCA as a preprocessing step, although we did estimate more general networks and found they performed no better.This could be a consequence of the data which were entirely observational, obfuscating causality.Understanding interventions using our results is similarly subject to the caveat that biomarkers behave the same whether they are observed or intervened upon 34 .This seems plausible since observational studies include everyday interventions such as disease, medicine, or lifestyle changes.Finally, our model is at the population-level and hence we cannot resolve homeostatic changes at the individual level.
We see exciting opportunities for future work.Our observation that principal components could be effectively used as independent variables suggests that more complex statistical models could also be applied.For example, individual-level model parameter estimates via mixed-effects modelling would help to determine whether individual health changes are gradual or sudden (or possibly critical 35 ).Changes in our model parameters due to age, chronic or acute illness, or medical interventions is particularly interesting, but will require specialized datasets to assess.Fortunately, small datasets are tractable with our linearized model.The generic nature of the model and its ability to find accumulating natural variables could also be applied to other biological or temporal scales.Others have postulated that damage aggregates due to dysfunction in regulatory systems or other intermediate scales 3 , which could be tested.Composite health measures, including biological age 31 , are also interesting to explore using our approach.Applying our approach using multiple biological ages as biomarkers 5 will naturally extract salient information regarding stability and mallostasis, as well as a smaller set of essential natural variables.New datasets will open up new opportunities for this analysis pipeline.It is interesting to consider leveraging the effect of the natural variables to intervene and observe in clever ways.For example, z 1 appears to be a biomarker of frailty, which affects both mental and physical health 32 , hence we could potentially intervene based on a physical mechanism but monitor using cognitive changes.
We have developed and applied a lightweight network model that includes the salient features of homeostasis: equilibrium values and recovery rates.Equilibrium values are allowed to drift, to accommodate allostatic changes.Across datasets and species we consistently observed that the linear decline of biomarkers with age was governed by a small set of accumulating natural aging variables.This accumulation can be described as mallostasis: homeostatic dysfunction and associated declining health.These variables appear to be important measures of age-related decline, including health and mortality.Their effects are spread out by a network of interactions, driving drift in the observed biomarkers, and potentially diluting and obfuscating the effects of age.We find that generic biomarkers spontaneously move towards an equilibrium position which is itself continuously drifting towards ill-health.Mallostasis is a generic feature of the aging process.
The Paquid dataset is a random subset of 500 humans (212 males and 288 females) from the Paquid prospective cohort study, enriched in dementia prevalence 36 .Age range: 66-95 years-old.Individuals were measured on average every 3.2 years for a maximum of 9 timepoints.We modelled four ordinal variables, including three measures of mental acuity: mini-mental state examination (MMSE), Benton visual retention test (BVRT) and Isaacs set test (IST), along with a self-reported depression score (CESD).We considered for covariates: sex, age and education level (completed primary vs not).
The Study of Longitudinal Aging in Mice (SLAM) includes two datasets, one for each mouse strain.Both include body composition measures and glucose serum at 12 week intervals starting at 7 weeks of age and continuing for the lifespan of each mouse 37 .Body composition and serum measurements were staggered and had to be imputed.Covariates included age and sex.We dropped 538/66138 measurements that were recorded after death, ostensibly these were coding errors.After preprocessing, the first dataset included 608 C57BL/6 mice (303 male and 305 female) measured on average every 6.2 weeks for a maximum of 20 timepoints (every 4.9 human-equivalent years).C57BL/6 mice are genetically similar (inbred) and prone to lymphoma and metabolic dysfunction 39 .The second included 611 Het3 mice (304 male and 307 female) measured on average every 4.2 weeks for a maximum of 27 timepoints (every 3.6 human-equivalent years).Het3 mice are a genetically heterogeneous 9/43 cross of four inbred mice (including C57BL/6) 39 .We converted to human-equivalent years using the ratio of median survival times of each strain to ELSA.Full details of the study are described elsewhere 24,37 .
The English Longitudinal Study of Ageing (ELSA) is a representative sample of English people aged 50 and over (with some younger) 38 .We used physical functioning questionnaire data and blood tests for 9330 humans (4063 males and 5267 females), reported at 4 timepoints, each separated by approximately 4 years.Our choice of 25 variables includes frailty measures, cardiometabolic biomarkers, and immune biomarkers (Supplemental Table S1).We considered waves 2, 4, 6 and 8, since only these contained the full suite of biomarkers.Covariates included age and sex.We considered only individuals whom were present both in wave 2 and in subsequent waves, thus excluding new recruits.Despite the large number of individuals, ELSA appeared to have the worst quality data due to high individual heterogeneity and low number of timepoints.

Data handling
All missing data were imputed.Dead individuals were also imputed, as it reduced bias due to mortality in simulated data (Supplemental Section S4.1).We compared several imputation strategies, including carry forward/back, multivariate imputation using chained equations (MICE) 40 , and using our model to impute the model mean.Ultimately, we used carry forward/back followed by the model mean, except for ELSA which used each individual's mean biomarker value followed by the population multivariate normal mean then model mean.See Supplemental Section S4 for details.
Estimation of our model, equation ( 1), used Supplemental Algorithm S1, which iteratively (×5) applied the maximum likelihood estimator: for Λ Λ Λ (which includes µ 0 through x 0 = 1), and for W W W , where the expectation values are to be taken over times, n, and individuals, i.For the diagonal models we instead used weighted linear regression.Missing values were imputed with the model prediction after each iteration (except ELSA).
Estimators are described and validated in Supplemental Sections S6 and S7, respectively.We used a time-dependent Cox model to assess survival.We assumed stepwise constant covariates via start-stop formatting 41 .All correlations are Pearson.All errorbars are standard errors unless stated otherwise.

Model assessment
We simultaneously estimated both parameter uncertainty and model performance using the standard deviation from 100× repeat bootstrap resampling.We compared model performance using the root-mean squared error (RMSE) and mean absolute error (MAE).Both were estimated using out-of-sample bootstrap 42 .In validation tests we found that a simple 632 estimator i.e., RMSE 632 ≡ 0.632 • RMSE test + 0.368 • RMSE train , provided a good estimate for the true values of both performance metrics (Supplemental Figure S7).0.632 is the expected fraction of unique individuals in each bootstrap 42 .

S1 Supplemental Information
We modelled generic health biomarker data as a mean-reverting stochastic process.Our study pipeline is summarized in Figure S1.Our model describes generic dynamics near an equilibrium solution (Section S8.4).In this supplemental we provide additional information to support and validate both our methods and our conclusions.We provide a complete description of the data in Section S2 and how we preprocessed it in Section S3.We performed a number of consistency checks on missing data which were imputed according to Section S4.We consider variations of our model in Section S5 which demonstrates that our final model best describes the data.We provide the mathematics necessary to estimate model parameters, along with an iterative estimation algorithm in Section S6.We then validate our algorithm using synthetic data in Section S7.Additional mathematics useful for understanding our model and its connection to the literature are described in Section S8.Finally, we include additional results in Section S9 which support our conclusions.• Scale by baseline standard deviation.

Figure S1
. Study pipeline.We analysed four datasets using our proposed model.We model the dynamics of biomarkers, ⃗ y n , over time using equation (S4).Our model extracts an interaction network, W W W , and equilibrium positions ⃗ µ n , where the latter are allowed to depend on covariates (e.g.age and sex).The estimated network, W W W , captures arbitrary linear interactions between biomarkers which can be removed by working with the natural variables,⃗ z n .Natural variables are defined by a linear mapping into the eigenspace of W W W .The natural variables allowed us to analyse stability.We were also able to infer changes to the mean and variance of the observed variables based on changes in the natural variables.We analysed 4 datasets derived from 3 longitudinal studies.The datasets and predictors ("biomarkers") used are summarized in Table S1.All predictor variables were continuous or, in the case of Paquid, ordinal with many scales (> 15).We included covariates to reduce confounding effects and to look for allostasis, which depends on age.We included age (continuous) and binary variables.The covariates used are summarized in Table S2.

S3 Preprocessing
The data we analyzed were longitudinal with regular sampling rates.For this reason, data were conveniently stored as 3-dimensional arrays (individuals, biomarkers, time points), meaning that each individual had the same number of variables and measurements (although many of them missing).This means that some timepoints for some individuals had to be 'invented' (instantiated as NA) based on the sampling rate of the study in question.
The Study of Longitudinal Aging in Mice (SLAM) datasets were both processed using the same criteria.The initial data were downloaded and processed using the analysis script of another publication 24 .We then applied additional preprocessing as follows.Biomarkers were visually investigated for normality and deemed adequate.The sex-specific mean and standard deviation of the first measurement of each biomarker was used to center and scale all timepoints.Mice with less than 2 timepoints were excluded from analysis (about 1% of mice).Any observations made past the reported death age of each mouse were excluded from analysis (about 1% of observations).We excluded the first two timepoints from analysis because after encoding we found that approximately half of individuals had not yet had a body composition measurement (imputed values looked unrealistic).Missing timepoints -which occurred due to staggered data collection -were instantiated using a piecewise linear model between known observations.Data from SLAM and the other datasets were stored in 3-dimensional arrays, with missing values imputed according to Section S4.The final arrays were size: (608, 6, 22) for SLAM C57/BL6, and (611, 6, 29) for SLAM Het3 (individuals, biomarkers, time points).
The Paquid dataset we used is available as part of a software package 36 .Biomarkers were visually investigated for normality.To improve normality we transformed CESD by the square-root and MMSE by − √ 30 − MMSE where 30 is the maximum allowed score for the MMSE.All biomarkers were centered and scaled by their respective mean and standard deviation from the first timepoint.Missing timepoints were instantiated using a piecewise linear model between known observations.The final array was size (500, 4, 9); (individuals, biomarkers, time points).
The English Longitudinal Study of Ageing (ELSA) dataset is available from the UK data service (https://ukdataservice. ac.uk/).We analysed all of the waves which included lab work: 2, 4, 6 and 8 (i.e. the "nurse" waves).We included only individuals present in wave 2, thus excluding later recruits.Biomarkers were visually investigated for normality.We found that the log transformation improved normality for C-reactive protein, ferritin, triglycerides, and white blood-cell count.All biomarkers were centered and scaled by their respective mean and standard deviation from the first timepoint.Skipped timepoints were instantiated using linear interpolation of the available timepoints.Censored (or died) timepoints were instantiated using the mean followup time (which was uniform due to the study design).The final array was size (9330, 23, 4); (individuals, biomarkers, time points).

S4 Missing data
We were presented with two forms of missing data for an individual at a particular timepoint.The entire timepoint could be missing or some subset of values could be missing.In either case the missingness could be informative; for example an individual may have temporarily left the study due to poor health and their biomarkers could have had abnormal values 13/43 reflecting their poor health.In this way the population may appear abnormally healthy as it ages.Under such circumstances, failure to impute can lead to biased study conclusions 45 , such as parameter estimates (Section S4.1).
We considered three imputation approaches and selected the approach which gave the most reasonable values.First ("simplest"), we imputed a single value using either the individual's mean biomarker value, carry forward the previous value (and then carry back any skipped values), the conditional population mean (assuming multivariate Gaussian statistics), or the individual mean followed by the conditional mean for individuals whom did not have that variable reported.Second ("model mean"), we considered an iterative approach after applying one of the first methods wherein values were imputed according to the model mean (i.e.model prediction), equation (S1) and equation (S3).Third ("MICE"), we considered multivariate imputation using chained equations (MICE) 40 .MICE is a multiple imputation technique that uses a Gibbs' sampler along with a predictive model.We considered MICE using both classification and regression trees (CART) and 2-level modelling (normal for continuous variables and logistic for binary).
When imputing the model mean, at each iteration we estimated the model parameters then imputed the conditional mean for each missing value (Algorithm S1).Let ⃗ y denote the biomarker, ⃗ u denote all unobserved y and ⃗ o denote all observed y.The statistics are Gaussian (equation (S4)), so we can use the factorization theorem 46 to compute the expectation value.
If ⃗ y n+1 is known but a set of ⃗ y n are unknown then where E(x|y) denotes expectation value of x conditional on y, Q Q Q is the precision matrix defined below, and I I I is the identity matrix.Note that is the block-decomposition of the noise covariance rearranged for observed (o) and unobserved (u) variables.If instead ⃗ y n is known but a set of ⃗ y n+1 are unknown then We used the simplest approach as an initial imputation (e.g.carry forward/back).We then imputed ⃗ y 1 first using equation (S1) then each subsequent timepoint using equation (S3).We compared the imputation quality and found that the model mean was both straightforward and effective, and therefore elected to use it for both SLAM datasets and Paquid.We initialized imputation with carry forward/back: that is, forward carrying previous values until the last timepoint was reached then backwards carrying to fill any values still missing.In rare cases a few data points were missing after imputation, these were simply ignored (we used all available case data).The ELSA dataset was sensitive to the model mean -perhaps due to the limited number of data points -hence we used a single imputation which combined first imputing the individual-specific variable mean followed by the conditional mean, assuming multivariate Gaussian statistics at each timepoint.Note that since we elected to use bootstrapping, we imputed each bootstrap replicate and then averaged to get an estimate for each missing value along with a standard error.
The final imputation was assessed for quality, Figure S2 -and looked reasonable.When inspecting imputation quality we are looking for the same age-dependent pattern for both the imputed and observed values, both in terms of mean and dispersion.Informative censorship is possible, so for variables with survival effects we can expect that missing values should be at higher risk because they include individuals whom were censored due to poor health (or death).Risk can be inferred by the direction of drift with respect to age: data points which look 'older' are likely higher risk.Hence imputed values may look a little 'older' than observed values.

S4.1 Informative censorship
Is it better to impute dropped individuals (dead, censored) or not?Dropped individuals may have abnormal biomarker values leading to their exclusion, i.e. informative censorship.There is "substantial" evidence that dropped individuals in longitudinal studies have worse health 48 and their health biomarkers will reflect this, leading to a potential survivorship bias.We used simulated data to test for potential bias and observed that -if done well -imputing values for dropped individuals can reduce this bias.
We simulated data from our model equation (S4) using randomly generated parameters then imposed informative censorship.We simulated 100 times with 100 individuals in each simulation.Each simulation included 2 biomarkers.Parameters and biomarkers were draw from normal random variables.The diagonal of W W W was mean −1/4, the off-diagonals were mean 0 and the overall standard deviation was 0.1.The mean µ 0 was 0 and standard deviation was 0.1.No covariates were simulated.The noise was diagonal, Σ Σ Σ = 0.5I I I (also used for instantiating the population).We censored using Gompertz statistics with proportional hazards for biomarker values 49 (shape: α = 0.1, scale: λ = 10 −5 ).The proportional hazards coefficients were randomly sampled from a normal distribution with mean 1 and standard deviation 0.1, this ensured that large values of the biomarkers were preferentially censored.
The results of the simulation are shown in Figure S3 for various imputation strategies, with the horizontal dashed line indicating unbiased results.We observed that a significant bias existed in both the diagonal and off-diagonal elements of W W W , which were systematically over-estimated if the data were not imputed.Conversely, if we used only the simple carry forward/back imputation, an even worse bias ensued in the opposing direction.If we used the model mean imputation, however, we reduced the bias in W W W to nearly 0 without significantly increasing bias in the other parameters.For this reason, we imputed all dropped individuals: censored and dead.

S6.2 Maximum likelihood estimators (MLEs)
We derive the MLEs for equation (S4) in full generality.For convenience define where we have p variables, N individuals and T + 1 timepoints.We index the N individuals with i and the T timepoint-pairs with n.For convenience we drop µ 0 and define the equivalent ⃗ µ in ≡ Λ Λ Λ⃗ x in ; where we use x i0n ≡ 1 to recover µ 0 .Estimators are denoted with a hat e.g.Ŵ W W estimates W W W .The log-likelihood is, We derive analytical forms for the MLEs as well as providing derivatives for gradient-based optimization algorithms.We also report the curvature since this is used to estimate the asymptotic error via the inverse Fisher matrix 51 .We found that the asymptotic errors are well-calibrated for W W W , but tend to be too small for ⃗ µ n (Section S7).In the present study, we report bootstrap errors.
Note that we find it useful to express the estimators in terms of the uncentered (cross)covariance, where the expectation value is taken over individuals, i, and timepoints, n.In general, ⟨ f (x in )⟩ i,n denotes the average of f (x in ) over individuals, i, and timepoints n.We start by considering W W W .The derivatives are where ∇ W denotes the gradient with respect to (vectorized) vec(W W W ). The MLE is thus The latter equation is useful for linear algebra software packages.Alternatively, we can invert the uncentered covariance the Fisher information is the negative of this.The covariance of the MLE is given by the inverse Fisher information, the standard errors are the square-roots of the diagonal elements,

S7.1 Parameter error
Here we test the calibration of our parameter errorbars.We compare both the bootstrap and asymptotic error estimates to the ground truth.Bootstrap errors were estimated using the standard deviation of bootstrap replicates.Asymptotic errors were estimated using the estimators in Section S6.2.As we will demonstrate in this section, the asymptotic errorbars can be too small, whereas the bootstrap errors appeared to be correctly calibrated.For this reason, we always used the bootstrap estimates in the main text.The asymptotic error estimates are much faster to compute and are presented for posterity.
In Figure S6 we present the coverage of both error estimators.The coverage is the fraction of times that the true parameter value fell within the estimated error interval.The nominal coverage is 68.3% for a normal random variable.In Figure S6A we present the coverage of the asymptotic error estimates and find that they are unsatisfactory for µ age and µ 0 (errorbars were too small).These may be due to strong correlations between the two parameters, for example the parameters for body weight correlated across simulations at cor(µ age , µ 0 ) = −0.886,which could make the asymptotic errors inaccurate.In Figure S6B we observe that all of the bootstrap error parameter coverages were close to the nominal rate (dashed line), and were symmetrically distributed above and below.This indicates that our bootstrap parameter errorbars were properly calibrated.

S7.2 Prediction error
Our primary measure of prediction error was the root-mean-squared error (RMSE).It is important that our measure is properly calibrated such that it estimates the correct RMSE, i.e. in a simulation study where the true error is known.We compare three RMSE estimators to the ground truth: (i) the testing error, which is the out-of-sample bootstrap error, (ii) the training error, which is the in-sample bootstrap error, and (iii) the 632 error which is a linear combination of 63.2% testing error and 36.8%training error 42 .The ground truth error is the error of the sample given the correct parameter values: this is the error of a single sample, not the distribution of possible values.The average ground truth error should be an unbiased estimate of the true, distribution error.In Figure S7A we demonstrate that the 632 error is close to the ground truth error.In Figure S7B we present the coverage of each estimator and find they are all close to the nominal rate.We conclude that the 632 error is a satisfactory estimator of the true model error.), the estimated value is plotted against the ground truth value for a variety of sample sizes (indicated by the legend).Points show mean; bands are the interquartile range (25th to 75th percentile).Bias is indicated by position of point relative to the red dashed line, y = x (perfect estimator).Precision (and accuracy) are inferred by the dispersion (bands).As the number of individuals, N, is increased from 50 to 1000 we see the estimator becomes increasingly accurate and precise, with a small dispersion around the ground truth values for each parameter.Points are staggered for visualization.Note: N = 10 had large errors and hence was excluded for better visualization.The coverage of the train and 632 error were close to the nominal rate, 68.3% (dashed line).The test error clearly had abnormally high coverage, indicating the errorbars on the test error are too large.Note: the true (stochastic) error is difficult to precisely estimate due to non-uniform sampling, so we used the average ground truth to estimate the true error.Errorbars are standard error.

S9 Additional Results
We restricted the main text to only our key results.Here we provide additional information to support our conclusions.We included covariates, ⃗ x, in the equilibrium position, µ(⃗ x), to reduce confounding effects and to test for the presence of allostasis (which depends on age).Here we tested each parameter for significance using the bootstrap parameter error estimates.The z-score for each covariate is reported in Figure S8; blue tiles are not significant, white and red are significant at p ≤ 0.05.Most covariates were significant, particularly age for the human studies.In Section S5 we found that the effect of covariates on prediction was typically small.This means that the effects of covariates were reliably estimated (small p) but did not explain much variation (minor effect on RMSE).
Our model estimates an interaction network, W W W , together with equilibrium positions.In the main text we presented the ELSA network with suppressed diagonal (Figure 2).The complete networks for each dataset are provided in Figure S9.The networks are all symmetrical because we used PCA as a preprocessing step.Relationships indicate how the y-axis variable will affect the x-axis variable during the next timestep.
Our model also estimates an equilibrium homeostatic position for each variable, µ.An important question is how strongly do variables adhere to homeostasis in the biomarkers, ⃗ y versus the natural variables, ⃗ z?In the main text we presented the difference between the natural variable mean and the equilibrium position for each variable, ⟨z j − µ j ⟩.We reproduce that figure beside the observed biomarkers in Figure S10.In Figure S10B (and Figure 3A) we observed that the natural variables appear to be split into two groups: the majority group was close to µ, indicative of homeostasis, whereas the minority group was far from homeostasis.This latter group had a strong drift term, µ age , which indicated that homeostasis was a moving target i.e. allostasis.In Figure S10A we show that the observed biomarkers were much more likely to be far from homeostasis than the natural variables (B), implying that the natural variables are able to condense the effects of age-related drift (allostasis) into a few variables (see also Section S8.6).
The natural variables appear to be efficient for representing age-related changes.What do the natural variables mean in terms of observable outcomes?In Figure S11 we report the correlations between the accumulating/drifting natural variables and biomarkers.In Figure S12 we report correlations with covariates.Together these give us an idea of what each natural variable represents and, by model implication, is controlling.For example, z 1 of Paquid is strongly correlated with the mental acuity scores: MMSE, BVRT and IST, implying it represents overall mental acuity.This may explain why z 1 was such a strong predictor of dementia (Figure S14).In this way the age-related decline of mental acuity can be represented by changes to just one variable, z 1 , but also observed across several biomarkers, MMSE, BVRT and IST.
The linear map, P P P, allows a few natural variables to cause several biomarkers to drift.The effects of allostasis on observed biomarkers via the primary risk natural variables are illustrated in Figure S13.The sign of each natural variable is arbitrary, due to idiosyncrasies in eigendecompositions 23 .The dominant survival dimension for the Het3 mouse data was z 2 , which appears to capture a loss of body fat and muscle, and relative gain of fluid.The dominant z 2 dimension for the C57BL/6 was more specific to loss of fat (the z 1 signal for C57BL/6 was very similar to the z 2 signal for Het3, Figure S11).z 1 was the dominant dementia-free-survival dimension for Paquid, and captured a system-wide drop in mental acuity scores (MMSE, BVRT and IST), which likely captures cognitive decline associated with dementia.z 1 for ELSA appears to be related to frailty 23 , having its effects spread across many variables, especially those related to disability (eye, hear, FI ADL and FI IADL), physical condition (grip strength and gait speed), and self-reported health (SRH); note that higher is better for physical condition variables and worse for the other variables (eye, hear, SRH, FI, etc).In all cases the effects are strongest in the natural variables, which is ensured by the orthogonality of P P P.This means that the effects of the natural variable drift must be diluted across the observed biomarkers (e.g. Figure S10), potentially hiding them within healthy variation.
The drift rate of the natural variables, µ age , was correlated with the risk of adverse outcome (Figure 4A).We named this phenomenon "mallostasis": the tendency of an aging system towards an ever-worsening equilibrium.Here we consider the role of confounding variables by constructing complete survival models for each natural variable and adverse outcome (mortality or dementia onset).We constructed a (time-dependent) survival model for each natural variable, with age, sex and the natural variable as predictors.We then recorded the Cox proportional hazards coefficients, which represent the (conditional) log-hazard ratios per unit increase for each natural variable.We observed that the Cox coefficients correlated with the drift rate, µ age : Figure S14.This provides more robust support for mallostasis: that the steady-state behaviour of aging mice and humans is declining natural variables and commensurately declining health.
As an illustration of mallostasis, we consider a simple composite health measure, b ≡ ⃗ µ T age ⃗ z. Figure S14 demonstrates that Cox coefficient is proportional to µ age therefore b is proportional to the hazard.This is confirmed in Figure S15 which demonstrates that b for each dataset is a good predictor of survival (or dementia onset).
The natural variables are connected to survival via mallostasis, but how do they relate to the observed phenotype?That is, how do the changes in the natural variables with age affect the observed biomarkers?The mean and variance are conserved between the biomarkers and the natural variables by Parseval's theorem.This means that natural variables with large means and variances will dominate the means and variances of the observed biomarkers, thus controlling the major changes we see.The

29/43
steady-state mean grows indefinitely proportional to µ age , what about the variance?The equilibrium (steady-state) dispersion, equation (S48), for the natural variables are plotted in Figure S16.Smaller eigenvalues are associated with higher variance.Some dimensions (e.g.z 1 for the C57BL/6) can contain as much as 10x more variance than the next highest dimension.These dimensions will dominate the observed variance in the steady-state.The model predicts that these dimensions will eventually become the dominant principal components (PCs), equation (S49), implying they would dominate the observed phenotype in the steady-state.Hence what we observe will be dominated by natural variables with small λ and large µ age , such as z 1 of ELSA, which appears to be closely related to frailty.
We used PCA (principal component analysis) as a preprocessing step, which allowed us to fit a diagonal model, equation S6.This simplified analysis and yielded equivalent performance to the full model (Section S5).Here we test the self-consistency of the approach.By using PCA, at each bootstrap the eigenvectors of W W W are principal components, possibly reordered (because we fit a diagonal model for W W W ). Averaging over multiple bootstrap replicates removes this equivalence -although in the steady-state the model predicts that the principal components and eigenvectors of W W W will coincide, equation (S49).Here we test the similarity of the PCA rotation and the eigenvector rotation: if they coincide then the principal components are eigenvectors.In Figure S17 we present the inner product between these matrices, which varies from −1 to 1, with ±1 representing perfect similarity.We observed strong similarities between the transformations, indicating that the principal components and natural variables will be strongly correlated.This suggests that PCA may be a useful shortcut for approximating the eigenvectors of W W W .
Finally, we include a survival summary for each dimension in terms of conditional Cox regression and the C-index in Figure S18.The values are identical to those used in Figures S14 (Cox coefficient) and 4 (C-index).This permits the reader to investigate the relative importance of each dimension.Comparing to the correlates of each dimension, Figure S11, one can infer potential mechanisms.For example, z 2 of C57BL/6 has a strong survival effect (low is bad) and shows increasing glucose and fat, which could indicate metabolic dysfunction, which C57BL/6 are prone to 39 .  .Tile colour indicates interaction strength (saturation) and direction (colour) of the interaction from the y-axis variable to the x-axis variable.Inner colour indicates the limit of 68% confidence interval (CI) closest to zero (i.e. standard error).Non-significant interactions, at 68%, have been whited-out.Variables are sorted by diagonal strength (increasing rate).The matrices are real and symmetric because the data were diagonalized by an orthogonal matrix (PCA).Figure S10.Homeostasis of biomarkers vs natural variables.The dysruption of homeostasis seems to be diffuse across biomarkers whereas it is concentrated into a few natural variables.A. Observed biomarkers were typically far from equilibrium (dotted line).B. In contrast, most natural variables were close to equilibrium.We inferred that variables close to equilibrium were in homeostasis whereas those far from equilibrium were allostatic.Together these plots suggest that the natural variables were able to condense the effects of allostasis into a few major variables. .This provides further information about what information each natural variable, z, contains.We expect the strongly drifting variables to exhibit correlations with age, though the sign of each z is arbitrary.Male is a binary sex indicator (1: male, 0: female); sex is the converse (0: male, 1: female).CEP is educational attainment level (1: attained primary, 0: did not).We consider the drift of the primary risk natural variables: z 1 for ELSA and Paquid and z 2 for SLAM.We observe a continuous drift in the natural variables.We also plot the drift of the biomarkers which is directly caused by each z via P P P. In this manner, a few natural variables can drive drift across several biomarkers.Since P P P is orthogonal (length-preserving) the drift of each natural variable must be diluted across biomarkers (at most a single biomarker can drift at the same rate).See also the correlation matrices, Figures S11 and S12.For the SLAM datasets we've included only timepoints where the average age was over 80 weeks.
. Allostasis drifts towards the risk direction.We fit a Cox model for each natural variable including age and sex as covariates.The Cox coefficient -i.e.log-hazard ratio (HR) per unit increase -correlates with the steady-state drift rate, µ age .The dominant risk direction for each dataset has been labelled by eigenvalue rank (e.g.z 1 is 01).The equilibrium standard deviation provides a native scale for each variable. .Shown are the dot products between the principal component rotation and P P P. The dot product assesses similarity between the transformations ranging from 1: identical, 0: orthogonal, and −1: identical with opposing sign.Identical transformations will generate identical natural variables.If the transformations are identical then all values on the diagonal should be ±1 (sign is arbitrary 23 ).We see that the dot products are often close to ±1, indicating that the transformations are very close, although they do not perfectly coincide.For each dataset, the top row corresponds to the Cox coefficient standardized by the equilibrium dispersion (ln (HR)/SD(z) e ) while the bottom is the C-index centered to 0 (C − 0.5).A Cox coefficient greater than 0 indicates that higher values are at increase risk and vice versa.A centered C-index greater than 0 indicates that higher values are at reduced risk and vice versa (opposite of the Cox coefficient).The Cox model is conditioned on age and sex (the same as Figure S14); the C-index is unconditioned.We see that in humans, the first dimension is the dominant determinant in risk of death (ELSA) or dementia (Paquid).It is less clear in mice, where allostatic drift is a better way to identify important survival dimensions (Figure S14 or Figure 4A).Inner colour indicates the limit of 95% confidence interval (CI) closest to zero (non-significant are red on blue or blue on red).

Figure 2 .
Figure2. A. ELSA interaction network.Tile colour indicates interaction strength (saturation) and direction (colour) of the interaction from the y-axis variable to the x-axis variable.Inner dot colour indicates the limit of the 95% confidence interval (CI) closest to zero (more visible point indicates lower significance).Non-significant interactions have been whited-out.Diagonal has been suppressed for visualization (see dotted lines in B).The matrix is real and symmetric because the data were diagonalized by an orthogonal matrix (PCA).Variables are sorted by diagonal strength in both A. and B. (increasing rate).B. Recovery rates in human-equivalent (h.e.) years i.e. negative eigenvalues (−λ ).The smallest recovery rates determine system stability21 .A recovery rate of 0.025 implies 1 − e −1 = 63% recovery after −λ −1 = 40 years (95% recovery after 120 years).The survival data all have similar minimum rates near 0.025, whereas the dementia data was faster (Paquid).The dotted lines are network diagonals (−W j j ); the solid lines are rates (−λ j ).

1 <Figure 3 .
Figure3. A. Position relative to equilibrium vs recover rate.Most natural variables were homeostatic (near equilibrium at 0).Some (labeled) variables were observed to be far from equilibrium; variables are labelled by rank e.g.01 ≡ z 01 has the fastest recovery (furthest left).B. Characterization of natural variable deviations from equilibrium using equation(8).Observe that ELSA is the only dataset where memory may dominate the system behaviour (ratio ≲ 1 = 10 0 ), indicating that the followup period may have been too short to reach a steady-state.In both figures only mouse (SLAM) data points over age 80 weeks were used since biomarkers had a u-shaped curve over the lifespan24 .

•
Center by baseline mean.

Figure S2 .
Figure S2.Final imputation quality check, visualized using principal component 1. A. C57BL/6 mice (SLAM).B. Het3 mice (SLAM).C. Paquid (human, dementia).D. ELSA (human).Imputed values appear to be reasonable for each dataset.Principal component analysis (PCA) was applied to each dataset in the entirety, flatted across timepoints.Good quality imputation (blue triangles) should show the same trend and dispersion as the observed data (red points).Censored individuals likely have worse health, so imputed values may look a little 'older' than the observed.Age-dependence is indicated by the solid lines with confidence intervals (cubic spline; geom_smooth with defaults47 ).Outlying points are highlighted (l ± 3 where l is the ordinary linear regression model).Data points were labelled as imputed (blue triangles) if the preponderance of the rotation weights were missing: ∑ i=missing |U i1 |/(∑ j |U j1 |) > 0.5; where U U U is the PCA rotation matrix.

Figure S5 .
Figure S5.Algorithm S1 validation.For the indicated parameters in each measurement (A.-F.), the estimated value is plotted against the ground truth value for a variety of sample sizes (indicated by the legend).Points show mean; bands are the interquartile range (25th to 75th percentile).Bias is indicated by position of point relative to the red dashed line, y = x (perfect estimator).Precision (and accuracy) are inferred by the dispersion (bands).As the number of individuals, N, is increased from 50 to 1000 we see the estimator becomes increasingly accurate and precise, with a small dispersion around the ground truth values for each parameter.Points are staggered for visualization.Note: N = 10 had large errors and hence was excluded for better visualization.

Figure S7 .
Figure S7.Bootstrap error calibration.632 error is a satisfactory estimator of the true error.A. Test error (out-of-sample) was biased high, training error (in-sample) was biased low, whereas 632 error was nearly unbiased relative to the ground truth.B. The coverage of the train and 632 error were close to the nominal rate, 68.3% (dashed line).The test error clearly had abnormally high coverage, indicating the errorbars on the test error are too large.Note: the true (stochastic) error is difficult to precisely estimate due to non-uniform sampling, so we used the average ground truth to estimate the true error.Errorbars are standard error.

Figure S11 .Figure S12 .
Figure S11.Natural variable correlates -biomarkers (predictors).A. C57BL/6 mice (SLAM).B. Het3 mice (SLAM).C. Paquid (human, dementia).D. ELSA (human).This helps to describe what information is in each natural variable, z, and therefore what each natural variable is capable of controlling.The sign of each z is arbitrary due to idiosyncrasies of the eigendecomposition.

Figure S13 .
Figure S13.Natural variable drift drives biomarker drift.A. C57BL/6 mice (SLAM).B. Het3 mice (SLAM).C. Paquid (human, dementia).D. ELSA (human).We consider the drift of the primary risk natural variables: z 1 for ELSA and Paquid and z 2 for SLAM.We observe a continuous drift in the natural variables.We also plot the drift of the biomarkers which is directly caused by each z via P P P. In this manner, a few natural variables can drive drift across several biomarkers.Since P P P is orthogonal (length-preserving) the drift of each natural variable must be diluted across biomarkers (at most a single biomarker can drift at the same rate).See also the correlation matrices, FiguresS11 and S12.For the SLAM datasets we've included only timepoints where the average age was over 80 weeks.

Figure S17 .
Figure S17.Principal components are very similar to the natural variables.A. C57BL/6 mice (SLAM).B. Het3 mice (SLAM).C. Paquid (human, dementia).D. ELSA (human).Shown are the dot products between the principal component rotation and P P P. The dot product assesses similarity between the transformations ranging from 1: identical, 0: orthogonal, and −1: identical with opposing sign.Identical transformations will generate identical natural variables.If the transformations are identical then all values on the diagonal should be ±1 (sign is arbitrary23 ).We see that the dot products are often close to ±1, indicating that the transformations are very close, although they do not perfectly coincide.

Figure S18 .
Figure S18.Survival summary. A. C57BL/6 mice (SLAM).B. Het3 mice (SLAM).C. Paquid (human, dementia).D. ELSA (human).For each dataset, the top row corresponds to the Cox coefficient standardized by the equilibrium dispersion (ln (HR)/SD(z) e ) while the bottom is the C-index centered to 0 (C − 0.5).A Cox coefficient greater than 0 indicates that higher values are at increase risk and vice versa.A centered C-index greater than 0 indicates that higher values are at reduced risk and vice versa (opposite of the Cox coefficient).The Cox model is conditioned on age and sex (the same as FigureS14); the C-index is unconditioned.We see that in humans, the first dimension is the dominant determinant in risk of death (ELSA) or dementia (Paquid).It is less clear in mice, where allostatic drift is a better way to identify important survival dimensions (FigureS14or Figure4A).Inner colour indicates the limit of 95% confidence interval (CI) closest to zero (non-significant are red on blue or blue on red).

Table S2 .
Covariate Summary Model selection.A. C57BL/6 mice (SLAM).B. Het3 mice (SLAM).C. Paquid (human, dementia).D. ELSA (human).Lower error is better.y-axis is 632-RMSE on left and 632-MAE on right.Horizontal lines indicate the best performing model.We are looking for the simplest model that consistently hits those lines across datasets.We considered models significantly worse if they do not have an error interval overlapping this line; prioritizing RMSE.