Generating synthetic aging trajectories with a weighted network model using cross-sectional data

We develop a computational model of human aging that generates individual health trajectories with a set of observed health attributes. Our model consists of a network of interacting health attributes that stochastically damage with age to form health deficits, leading to eventual mortality. We train and test the model for two different cross-sectional observational aging studies that include simple binarized clinical indicators of health. In both studies, we find that cohorts of simulated individuals generated from the model resemble the observed cross-sectional data in both health characteristics and mortality. We can generate large numbers of synthetic individual aging trajectories with our weighted network model. Predicted average health trajectories and survival probabilities agree well with the observed data.


S1.2 Estimating likelihood
To calculate the likelihood in Eqn. [7], we run S simulations with the same set of parameters (with S > 10 6 individuals). We use a discrete kernel density estimate of the health term of the likelihood, where S(t (m) ) is the number of simulated individuals alive at age t (m) . The superscript (s) indicates a simulated individual and (m) an observed individual from the data. We use the Racine and Li kernel [4], otherwise.
This allows individuals that differ in damage by 1 deficit to contribute to the likelihood, as determined by an individual bandwidth λ (m) . The bandwidth λ (m) is selected by minimizing the mean-squared error [5], where p (m) is the empirical frequency estimate of the distribution (obtained using bandwidth λ = 0). The bandwidth λ (m) is individual-dependent, so that the well sampled individuals will have a small bandwidth, and poorly sampled individuals have a larger bandwidth to reduce noise. Both the censored and uncensored mortality terms of the likelihood are calculated by binning the death ages of simulated individuals that match the observed data ( i |d

S1.3 Regularization
To have an increase in mortality rate with decreasing health, the deficit mortality contributions β j in Eqn. [4] should all be positive. Without such a bound, we have observed that the optimization sometimes converges to positive and negative β j values scattered around zero, which leads to the deficit contribution to the mortality rate being negligible, and gives uniform mortality for all individuals. However in contrast, setting a strict bound β j ≥ 0 causes the optimization to converge to parameters which perform poorly with respect to the health trajectories. Instead, we have found that a soft penalty behaves well when added to the log-likelihood to penalize negative values, +C j min(0, β j ). (S4) Values of β j ∼ 1 give a significant mortality contribution, and log-likelihood values are estimated with a stochastic noise of ±10, so we choose C = 100 by hand to achieve a moderate regularizing effect -and we find that the optimization performs well with this choice.

S1.5 Parameter optimization
Optimization is done with particle swarm optimization (PSO) [6]. PSO is derivative-free and highly parallelizable. We use the standard version: The parameter values for the ith particle at an iteration t are represented as θ i,t and the velocity for this particle as v i,t . The current best set of parameters found by particle i is θ p i,t , and the current global best set of parameters is θ g t . We randomly sample u p , u g uniformly from [0, 1]. We set ω = 0.7, c p = c g = 2 and use 100 − 200 particles, depending on compute resources available. The optimization typically converges in 100-200 iterations.
Parameter bounds from section S1.4 are implemented by rescaling the velocity components v i,k,t with a hyperbolic method so that parameters never exceed their upper bounds U k and lower bounds L k [7], Our objective function is stochastic, so to avoid optimizing to the noise, every 10 PSO iterations we recalculate the global maximum likelihood for the global best set of parameters θ g t , and the individual particle maximum likelihood values for each particle's best set of parameters θ p i,t .

S2 Parameter robustness
To examine parameter robustness, we select new parameter sets by randomly selecting 25% of parameters to randomly perturb from their maximum likelihood valuesθ d within the range Broadly following the approach of [8], we weight these perturbed parameters by calculating an approximate log-likelihood by using only the χ 2 of first and second order joint distributions of health and population survival. This amounts to a Gaussian approximation to the log-likelihood, We use this estimate of L to efficiently estimate any f (θ): In Fig. S1 we show deficit prevalence and survival averaged over likelihood-weighted parameterizations via Eqns. S9 and S10. These average predictions are quite close to the maximum likelihood estimates, which indicate that prevalences and mortality are robust predictions of our modelling approach.
In contrast, simple measures of the network structure do not appear to be robust. In Fig. S2 we perform the optimization 13 times, and the compare network degrees for the 10 nodes of our network. We find a broad ranges of degrees for each node, with significant overlap of these ranges between nodes. This indicates that though the behavior of the networks is similar (see e.g. Fig. S3), the network structures themselves are not robustly predicted by the available data.

Preparing meal difficulty
6. House work difficulty 7. Take medicine difficulty 8. Managing money difficulty of node i given that node j is damaged, vs age. Each curve is averaged over combinations of the other node states, weighted by their probability of occurring in the simulation. Each intersecting row and column represents the indicated node labelled on the diagonal, and the direction of the links is indicated by labels on each subplot. Each subplot has 13 rate curves, each with different estimated parameters for a different starting seed. The qualitative agreement of these curves for particular pairs of nodes indicates that these detailed damage rates are qualitatively robust (note the linear scale), even though the effective degrees of the nodes are not (see Fig. S2). Observed NHANES prevalence is shown with red circles; standard errors are smaller than the point size. Here predictions do not always start at value from the observational data (comparing the agreement for the first red point with Fig. S4) because there is missing data that we impute by simulating from birth.   Distribution of the Frailty Index, F = N i=1 d i /N for different age ranges for the model (unfilled blue histograms) and the CSHA data (filled red histograms). Note that the maximum FI for small numbers of deficits (in this case N = 10) is 1, as is reached in the CSHA data at later ages. Recovering this FI limit with model data, as well as broader consistency with the CSHA distributions, confirms that we capture the heterogeneity of the population as individuals age.  (7,9) 70 80 90 (8,9)      We see that the correlation structure in all three sets of deficits is quite different. We see stronger positive correlations in the CSHA deficits used for the main plots than the alternative CSHA deficits. We generally see positive correlations for the NHANES deficits, but see a negative correlation with all deficits for the deficit "Coronary heart disease". This deficit also has slightly decreasing prevalence with age, so this is likely due to those individuals with coronary heart disease dying early. In D), we show the probability of an inferred weight in the model being positive, for 13 different model fits for the same deficits as A). The sign of the weights is quite different from the correlations. Probabilities near 0.5 mean that the inferred connection is not very robust.