Synchrony of biomarker variability indicates a critical transition: Application to mortality prediction in hemodialysis.

Critical transition theory suggests that complex systems should experience increased temporal variability just before abrupt state changes. We tested this hypothesis in 763 patients on long-term hemodialysis, using 11 biomarkers collected every two weeks and all-cause mortality as a proxy for critical transitions. We find that variability-measured by coefficients of variation (CVs)-increases before death for all 11 clinical biomarkers, and is strikingly synchronized across all biomarkers: the first axis of a principal component analysis on all CVs explains 49% of the variance. This axis then generates powerful predictions of mortality (HR95 = 9.7, p < 0.0001, where HR95 is a scale-invariant metric of hazard ratio; AUC up to 0.82) and starts to increase markedly ∼3 months prior to death. Our results provide an early warning sign of physiological collapse and, more broadly, a quantification of joint system dynamics that opens questions of how system modularity may break down before critical transitions.

Introduction

Apparently stable systems can sometimes undergo abrupt shifts, known as critical transitions. Critical transitions can reflect either a generalized collapse of the system, such as the collapse of the cod fisheries and associated ecosystems in the western North Atlantic in the 1990s or financial crises, such as occurred in the US in 2007–2008 (May et al., 2008). Early warning signs (EWSs) of critical transitions could allow for corrective action, particularly to avoid system collapse (Macrae, 2014). EWSs have been studied in a number of contexts, particularly in ecology (Scheffer et al., 2015), and may include increases in variance (Carpenter and Brock, 2006; Drake and Griffen, 2010), change in auto-correlation structure (Ghanavati et al., 2014; Gijzel et al., 2017), and increases in cross-correlation (Scheffer et al., 2009), all of which are part of a phenomenon known as “critical slowing down” that precedes the transition (Dakos and Bascompte, 2014). These indicators of change in system state are derived from structured changes in the dynamics at the onset of a critical transition.

Recently, there is a growing interest in applying critical transition theory to understanding human health and physiology in order to predict progression toward disease states (Almeida and Nabney, 2016; Chen et al., 2012; Gijzel et al., 2017; Trefois et al., 2015). In many medical contexts, apparently stable patients decompensate or decline relatively rapidly. Examples include epileptic seizures (Kramer et al., 2012), decompensation in intensive care (Ghalati et al., 2019), and clinical frailty in the elderly, which represents a point of accelerated decline (Fried et al., 2021; Gijzel et al., 2017). Several studies have started to successfully apply critical transition theory in a medical context (Almeida and Nabney, 2016; Chen et al., 2012; Gijzel et al., 2017; Maturana et al., 2020; Trefois et al., 2015; van de Leemput et al., 2014), but few have fully drawn on the multivariate nature of biological systems. Multivariate indices of EWSs might successfully incorporate not just the dynamics of individual biomarkers preceding transitions but also the joint dynamics of the numerous interacting biomarkers (Cohen, 2016), substantially increasing predictive power.

We propose that beyond specific diseases (Li et al., 2014; Maturana et al., 2020; Rockne et al., 2020; Tarazona et al., 2019), the broad health status of an individual can be assessed using critical transition theory, with general outcomes such as death often reflecting a critical transition, which we refer to as a “physiological collapse” of the organism. Successful prediction of impending physiological collapse could allow timely interventions to reduce mortality and/or start end-of-life planning. However, one challenge in a medical context is obtaining sufficiently detailed time series data to calculate meaningful EWSs. Clinical data are often subject to reporting biases (e.g., measures taken when patients are sick, selective reporting of diagnostic codes), and cohort studies are rarely if ever at a fine enough temporal scale. An exception is hemodialysis data in patients with end-stage renal disease: in most contexts, hemodialysis occurs 3x/week with blood draws approximately bi-weekly, continuing over several years, until kidney transplant or death (Nakazato et al., 2017). Accordingly, data on patients under hemodialysis present a time series of largely complete biomarker data generally unavailable in other medical contexts, and permit us to assess changes in variability in ways that are generally not possible in cohort data or other standard clinical contexts. Recently, albumin variability was shown to increase before death in a Japanese cohort of patients under hemodialysis (Nakazato et al., 2017), and a multivariate index of variability predicted frailty in the same cohort (Nakazato et al., 2020).

Here, using mortality as a proxy for physiological collapse in patients under hemodialysis, we predicted that the variability of biomarkers would provide a stronger prediction of impending mortality than the levels of the biomarkers. We likewise predicted that the general behavior of the system, as characterized by joint (multivariate) signals of the markers, would give stronger predictions than the levels or variability of any individual marker. Accordingly, we expected integrative measures of variability to be powerful EWSs. We assessed (1) whether there is a generalized joint signal of biomarker variability (i.e., synchrony of variability), and (2) if any such signal could predict critical transitions, in our case, death.

Results

Study population and biomarker selection

We studied biomarker trajectories in all patients with chronic kidney disease (CKD) on long-term dialysis in the Eastern Townships region of Quebec, based on data extracted from a clinical research database containing electronic patient records (see STAR Methods for details). After exclusions, our study population consisted of 763 patients (Table 1, Figure 1) followed for a total of 56,019 visits with biomarker data over 2,496 patient-years (median [interquartile range]: 45 visits/patient [15, 101] and 2.1 years/patient [0.7, 4.4]). Patients excluded from analyses were less likely to have diabetes, die, or have a kidney transplant, and were followed for a shorter time period (see Table S1). 525 patients died during follow-up, but 91 of these had no visit within 3 months prior to their death, likely due to transfer to palliative care. We detected 391/763 (51.2%) of patients as diabetic. We generated two lists of biomarkers based on their measurement frequency: one that includes 11 biomarkers generally measured every two weeks (2-week list), and one including 16 biomarkers measured every four months (4-month list, see Table 1). Main analyses were conducted with the 2-week list, while we used the 4-month list to validate our results at a different time scale and incorporate albumin, which was only measured every four months in our cohort.

Table 1

Study population characteristics at first visit included in analyses

Characteristic	(n = 763)	2-week list	4-month list
Age (years), mean ± SD	64.2 ± 15.8
Age range (years), min - max	16.2 - 94.6
Men, n (%)	479 (62.8)
Diabetics, n (%)	391 (51.2)
Death, n (%)	525 (68.8)
Kidney transplant (KT)
Patients who ceased HD after having a KT, n (%)a	93 (12.2)
Total number of patients with KT, n (%)	186 (24.4)
Follow-up length (years), median (Q25, Q75)	2.1 (0.7, 4.4)
Albumin (g/L), mean ± SD	36.1 ± 5.9		X
Calcium (mmol/L), mean ± SD	2.25 ± 0.19		X
Creatinine (μmol/L), mean ± SD	667 ± 303		X
Glucose (mmol/L), mean ± SDb	8.1 ± 5.0		X
Hematocrit (proportion of 1), mean ± SD	0.32 ± 0.05	X	X
Hemoglobin (g/L), mean ± SD	108 ± 15	X	X
MCH (pg), mean ± SD	31.3 ± 2.2	X	X
MCHC (g/L), mean ± SD	333 ± 11	X	X
MCV (fL), mean ± SD	94.0 ± 6.1	X	X
Phosphate (mmol/L), mean ± SD	1.67 ± 0.53		X
Platelet count (109/L), mean ± SDc	215 ± 80	X	X
Potassium (mmol/L), mean ± SD	4.78 ± 0.74	X	X
RBC count (1012/L), mean ± SD	3.5 ± 0.5	X	X
RDW (%), mean ± SDb	15.5 ± 1.9	X	X
Sodium (mmol/L), mean ± SD	138 ± 4	X	X
WBC count (109/L), mean ± SDb	7.8 ± 3.3	X	X

Abbreviations: KT, kidney transplant; MCH, mean corpuscular hemoglobin; MCHC, mean corpuscular hemoglobin concentration; MCV, mean corpuscular volume; Q25, 25th percentile; Q75, 75th percentile; RBC, red blood cells; RDW, red cell distribution width; SD, standard deviation; WBC, white blood cell.

Within 30 days from the kidney transplant.

Variable was log-transformed.

Variable was square root-transformed.

Figure 1

Study population flowchart diagram

Of the 2,565 patients having hemodialysis (HD) at the CHUS between January 1997 and December 2017, we excluded 1694 patients who had HD for less than six months and 73 patients who had both an irregular HD visit frequency and an acute or unspecified kidney failure (KF) diagnosis. Then, we selected visits with biomarker data at least three days apart, and excluded nine patients with incomplete biomarker data (including four without any biomarker profiles) and 26 patients with less than three visits in total. The final study population thus included 763 patients, among whom 413 died within 30 days from their last HD visit and 112 who died more than 30 days after their last HD visit. Within patients still alive at the end of the study period, 59 had a kidney transplant (KT) after which they stopped HD (within 30 days), 117 stopped HD for unknown reasons, and 62 were still on HD within 30 days from the end of the study period. See Tables S1 and S2 for demographic characteristics of patients excluded from analyses and of observations with incomplete data, respectively.

Prediction of mortality by levels and variability of individual biomarkers

First, we tested the hypothesis that individual biomarker variability would better predict mortality than individual biomarker levels. Variability of individual biomarkers was calculated as within-individual coefficients of variation (CVs) per year before death or censoring. To ensure comparability (since CVs, unlike biomarker levels, cannot be calculated at a specific timepoint), levels of biomarkers were calculated as within-individual means per year. We thus hereafter refer to information on levels as “means,” just as variability is referred to as “CVs.” We compared indices using Cox proportional hazards models to extract HR95, a scale-invariant hazard ratio for continuous variables that assesses variation in hazard across the range of the variable (Milot et al., 2014). Figure 2A shows mortality predictions for each biomarker, either using its yearly mean or CV, as well as p-values of proportional assumption tests for each coefficient. All models control for age (using a cubic spline), sex, diabetes diagnosis, and length of follow-up, clustering multiple observations per individual. Overall, CVs were better in predicting mortality compared to means: in all 11 biomarkers, higher variability (CVs) was associated with increased mortality, and nine had HR95 ≥ 4 (Figure 2A). In comparison, nine biomarkers predicted mortality significantly when using their levels, but only one had HR95 ≥ 4 (red cell distribution width [RDW], HR95 = 4.96, p < 0.001; see Figure 2A). For 10 of 11 biomarkers, the HR95 estimated from CVs was greater than that estimated from levels (Figure 2A).

Figure 2

Prediction of death by biomarker levels and variability indices

(A and B), HR95, together with 95% confidence intervals, are shown for the levels (means, red) and variability (CVs, blue) of each biomarker considered (A) and integrative multivariate indices (i.e. each principal component calculated on all biomarkers, (B). Levels of hemoglobin, hematocrit, MCH, MCHC, MCV, potassium, sodium, and RBC were inversed (1/x) to obtain HR95 above 1, for ease of representation.

(C) Receiver operating characteristic curves for a basic model including only demographic and control variables (black) and for models sequentially adding PC1 (green), CVPC1 (red), and all 10 other CVPCs (blue) are shown.

(D) Accuracy of mortality prediction for the first principal component of a PCA performed on means (PC1, blue) or CVs (CVPC1, red), or on either one controlling for the other PCs/CVPCs in the cox model (darker hues), by sequentially increasing the number of PCs/CVPCs added in the cox model. Cox proportional hazard models were performed with (dashed lines) or without (solid lines) including demographic control variables (age, sex, diabetes diagnosis, and length of follow-up).

See Figures S1, S2 and S6 for results of sensitivity analyses and Figure S3 for effects of combining means and CVs on mortality prediction. See also Figures S7 and S10 for results with the 4-month variable list and Figure S13 for the effect of the number of observations included in CV calculation. Abbreviations: AUC, area under the curve; MCH, mean corpuscular hemoglobin; MCHC, mean corpuscular hemoglobin concentration; MCV, mean corpuscular volume; RBC, red blood cells; RDW, red cell distribution width; WBC, white blood cells.

Integration of individual biomarkers into multivariate indices

To test the hypothesis that integrative measures would outperform univariate ones, we then combined information from individual biomarkers into integrative indices using principal component analysis (PCA), both for biomarker means and CVs, and extracted principal components (hereinafter referred to as “PCs” for means and “CVPCs” for CVs). The first principal component for biomarker variability (CVPC1) was by far the strongest predictor of mortality (HR95 = 9.68, p < 0.001, after controlling for age, sex, diabetes diagnosis, and length of follow-up; see Figure 2B), although several other indices also provided reasonably strong predictions (e.g., PC6 for biomarker mean levels, HR95 = 4.89). Results were similar when excluding 112 patients who died but were not followed in the last 30 days before their death, with CVPC1 even slightly more predictive (HR95 = 11.49, p < 0.001, see Figure S1), or when excluding highly correlated variables (hematocrit, red blood cell count, and mean corpuscular volume, see Figure S2A). Compared to other predictive indices, CVPC1 was also the most stable across population subsets, with HR95 > 5 in all subgroups tested (Figure 3 and Figure S2B). Prediction between men and women did not differ significantly, except for PC6 (HR95 = 3.5, 95% CI [2.5, 4.9] and HR95 = 7.6, 95% CI [5.1, 11.3], respectively for men and women).

Figure 3

Mortality prediction by PC1, PC2, PC6, CVPC1, and CVPC3 in different population subsets

(A–E) HR95, i.e. the hazard ratio of being in the 97.5th percentile relative to the 2.5th percentile of the index, together with 95% confidence intervals are shown for PC1 (A), PC2 (B), PC6 (C), CVPC1 (D), and CVPC3 (E) in different subsets of the study population. All models control for age using a cubic spline (with 3 degrees of freedom when performed on a specific age group and five otherwise) and length of follow-up, clustering multiple observations per individual. Models also control for sex and diabetes diagnosis, except when population is stratified using this variable. P-values of proportional assumption tests for the given coefficients are indicated.

See Figures S2 and S13 for results of sensitivity analyses using non-redundant biomarkers and at least five observations in CV calculation, respectively. See Figure S9 for results with the 4-month list.

Next, we assessed whether prediction surpassed the one obtained from using solely demographic and control variables (age, sex, diabetes diagnosis, and length of follow-up). To do so, we plotted receiving operating characteristic (ROC) curves of a 1) “basic” model using only demographic and control variables, and 2) of additional models sequentially adding PC1, CVPC1, and finally all 10 remaining CVPCs (Figure 2C). The obtained ROC curves show model improvement when adding CVPC1, while area under the curve (AUC) is nearly unchanged by the addition of PC1 or of the 10 remaining CVPCs. We also assessed whether prediction could be improved by combining information from our two types of integrative indices (means and CVs), or by sequentially combining all PCs or CVPCs one by one (Figure 2D). The AUC increased from 0.677 to 0.741 when adding all 11 PCs to CVPC1 in the model, and reached 0.813 when including covariates (sex, age, diabetes diagnosis, and length of follow-up). However, adding subsequent CVPCs to the model had only a minor effect on prediction accuracy (AUC = 0.822 compared to 0.813), with CVPC3 having the most effect (Figures 2C and S3). Broadly, the AUCs show that there are independent contributions of information from biomarker levels (means), variability (CVs), and demographic data, with the CVs contributing more than the means and CVPC1 providing by far most of the signal among CVPCs.

CVPC1 as a joint signal of biomarker variability

Because CVPC1 provided such strong prediction of mortality, we took a closer look at what it quantifies. CVPC1 captured an impressive proportion of overall system variability, explaining 48.8% of total CV variance, whereas the second and third principal components (CVPC2 and CVPC3) explained only ∼12% and ∼10% of the variance, respectively (Figures 4A and S2C). Variance was also more concentrated in the first axis for CVs relative to raw variables. Most strikingly, CVPC1 is driven jointly by all biomarkers rather than a subset: Loadings of individual biomarker CVs were nearly equal, a result highly replicable across men, women, diabetics, non-diabetics, and groups at various time intervals before death (Figure 4B). Supporting this, all pairwise correlations among biomarker CVs were positive, statistically significant at α = 0.05 (Figures 4C and S2D), and generally of similar magnitude (mean r = 0.35, SD = 0.18, min = 0.12, max = 0.96; see Figure 4D). These results show that variability is highly synchronized across the individual biomarkers: as one biomarker became more variable, all the others tend to be as well. In contrast, biomarker means were much less synchronized than CVs: PC1, which by definition explains the most variance of any axis in the analysis, explained only ∼25% of total variance in biomarker levels, and contributions of individual biomarkers were much less equally distributed (see Figures 4A–4C and S4).

Figure 4

Physiological variability shows a strong coordinated signal distributed evenly across all measured biomarkers

(A) Variance explained by PCA on raw biomarkers (triangles) or coefficients of variation (circles), for different population subsets relative to time of death or by demography.

(B) Relative biomarker contributions to CVPC1, ordered from largest contribution (hemoglobin) to smallest (MCHC) in the full dataset (see STAR Methods). Subsequent columns are based on loadings of the PCA run exclusively on the indicated subsets. Contribution for a given biomarker is the absolute value of the loading divided by the sum of the absolute values of all loadings.

(C) Pearson correlations (Corr) among raw biomarkers, coefficients of variation, and composite indices (CVPC1-3: First through third axes of the PC on coefficients of variation). Blue indicates positive correlations, and red represents negative correlations. Xs represent correlations not significant at α = 0.05. Above the diagonal are the CVs, and below are the biomarker levels.

(D) Histogram of pairwise Pearson correlation coefficients between CVs of individual biomarkers.

See Figure S2 for results of sensitivity analyses using non-redundant biomarkers and Figure S4 for variable contributions to PC1, PC2, PC6, and CVPC3. See Figures S8 and S11 for results with the 4-month variable list. Abbreviations: MCH, mean corpuscular hemoglobin; MCHC, mean corpuscular hemoglobin concentration; MCV, mean corpuscular volume; RBC, red blood cells; RDW, red cell distribution width; WBC, white blood cells.

Timing of change in CVPC1 before death

Finally, we assessed how long before death the joint variability increase could be detected by applying change point analysis to our integrative measure of biomarker variability (CVPC1) calculated with different time windows (i.e., every 2, 3, or 4 months). A stable change point was detected at ∼3 months before death, precisely at 3.3, 3.8, and 3.1 months, respectively, for CVPC1 calculated bimonthly, trimonthly, and four-monthly (Figure 5A). In comparison, change point analyses on other CVPCs yielded much larger confidence intervals, suggesting no clear change in their distribution before death (Figure 5B). This rise in CVPC1 before death, which is more abrupt than for any of the other PCs or CVPCs (Figures 5C and 5D), is also weakly discernible in randomly chosen individual profiles, and even prior to hospitalization events (Figure 6). Biomarker levels increase or decrease at different rates preceding death (Figure 5E), whereas increases in CVs show striking synchrony in the few months preceding death (Figure 5F). Among individual biomarkers, RDW levels show the greatest increase before death (Figure 5E). We had decided to exclude the first six months on hemodialysis based on the assumption that this period represents itself a critical transition (Broers et al., 2015). To validate this assumption, we looked at trends of indices after hemodialysis initiation by including all patients and observations, except patients with both irregular hemodialysis visits and an acute renal failure diagnosis (n = 1815 patients; see Figure S5). After excluding observations with missing data, we were left with 1765 patients with at least two observations. The general CVPC1 trend, calculated every three months, clearly shows an elevated overall variation in biomarkers in the year after hemodialysis initiation (Figure S5B), which is supported by trends of individual CVs (Figure S5D). Consequently, CVPC1 predictions of mortality are substantially lowered when considering this critical adaptation phase and including patients who died within it (Figure S6). Nevertheless, excluding only the first three months, instead of six, affects CVPC1 mortality prediction only moderately (HR95 = 8.2 compared to 9.7).

Figure 5

CVPC1 trend before death

(A) CVPC1, the first principal component of a PCA performed on CVs calculated every 2 (red), 3 (blue), or 4 (green) months, is plotted against time before death. Change point analyses were applied to regression models between CVPC1 and time before death, allowing slopes to vary across individuals, and are represented by the vertical dashed lines with the respective colors.

(B) Results from change point analyses applied to regression models between each CVPC (calculated every 3 months) and time before death are indicated with 95% confidence intervals.

(C and D) Integrative multivariate indices for biomarker levels (PCs, C) and variability (CVPCs, D) calculated every 3 months are averaged and plotted against time before death, centering at 5 years before death for ease of comparison. Results from change point analyses performed on CVPCs are indicated as vertical dashed lines.

(E and F) Biomarker levels (mean z-scores, E) and variability (CVs, F) were calculated every 3 months and averages are plotted against time before death, centering at 5 years before death.

See Figure S5 for index trends after hemodialysis initiation and Figure S12 for results with the 4-month list. Abbreviations: MCH, mean corpuscular hemoglobin; MCHC, mean corpuscular hemoglobin concentration; MCV, mean corpuscular volume; RBC, red blood cells; RDW, red cell distribution width; WBC, white blood cells.

Figure 6

Individual CVPC1 trajectories before death or censoring

The color represents the status at the end of follow-up (red for patients who died and blue for patients who were alive) and line type represents the diabetes diagnosis (a solid line for non-diabetics and a dashed line for diabetic subjects). Vertical green lines represent hospitalizations. Individuals were randomly chosen from among those with ≥6 years of follow-up.

Validation using the 4-month biomarker list

Albumin, whose variability was previously found to predict mortality in patients who were hemodialyzed (Nakazato et al., 2017), was only measured every four months in our study population. Therefore, we replicated key analyses using the 4-month list of 16 biomarkers (Table 1), for which a total of 11,475 visits were available, with a median [interquartile range] of 8 [3, 18] visits per patient. Results were broadly similar (Figures S7–S12), showing the stability of patterns regardless of the precise biomarker choice or temporal scale of the measurement. Results using the 4-month variable list showed smaller effect sizes for mortality prediction (HR95 = 4.99, p < 0.001 for CVPC1, Figure S7 and S9), but increase in variability among individual biomarkers still appears to be synchronized: CVPC1 explains ∼35% of the total variance in biomarker variability and is driven by all biomarkers rather than a specific subset (Figures S8A and S8B), and 117 out of 120 correlations among CVs are statistically significant at α = 0.05 and of roughly similar magnitude (Figures S8C and S8D). Albumin appears to show the greatest increase in variability in the year preceding death (Figure S12B).

Discussion

In this study, we have shown that biomarker variability predicts mortality during hemodialysis more strongly than biomarker levels, and that integrative indices are more predictive than individual biomarker variation. We also found a sudden, coordinated rise in the variability of all biomarkers, whose integration into a multivariate index is also a strong predictor of impending mortality: this index can detect 11-fold differences in mortality risk. Interestingly, the single biomarker showing the highest predictive performance and pre-mortem increase is RDW, a finding consistent with previous studies (Mucsi et al., 2014; Oh et al., 2012; Solak et al., 2014). RDW is itself a kind of CV—quantified through a simple equation (i.e. [standard deviation of red blood cell volumes] / [mean corpuscular volume] × 100)—and can thus be regarded as an indicator of momentary variability.

The coordinated rise in variability is striking in that even uncorrelated biomarkers from apparently unrelated systems show synchronized increases in variability (Figures 4 and 5F). Indirectly regulated (e.g. hemoglobin) and directly regulated (e.g. sodium) variables also appear to dysregulate in synchrony, and their increases in variation translate into mortality predictions of similar magnitudes (Figure 2A), contradicting other recent findings (Rector et al., 2021). In other words, the signal of the impending critical transition propagates quickly and diffusely across the whole system. While the biological underpinnings of this synchronization are not yet known, propagation of the signal across biomarkers is consistent with complex systems theory, and with recent findings suggesting that loss in resilience (frailty) emerges from parallel dysregulation in different physiological systems (Fried et al., 2021; Ghachem et al., 2021). Recent work has highlighted the need to understand whether critical transitions touch entire systems or only subsets of the systems (Weinans et al., 2021); our results appear to answer this question, at least in the case of hemodialysis, and for the measured biomarkers. More broadly, it will be important going forward to evaluate whether the synchronization of variability equates to a breakdown in modularity, and how general this is in critical transitions. The PCA approach used here is a promising method.

Our findings confirm that critical transition theory can apply to health, not only for highly specific physiological events such as epileptic seizure (Maturana et al., 2020) but also to overall health/homeostasis and its converse, physiological collapse. Here, we have shown applications of biomarker dynamics to patients with CKD under hemodialysis, but many other medical applications could readily be developed, including vital rates in intensive care (Ghalati et al., 2019), brain waves to predict sleep cycles (Bashan et al., 2012), and advent of clinical frailty (Gijzel et al., 2017). Also, with more regular data collection, there could be applications in congestive heart failure, predicting myocardial infarction, chronic obstructive pulmonary disease, and many others, as well as to fields such as ecology, economics, and climate science. Existing approaches to predict critical transitions/tipping points in health often rely on machine learning-based algorithms (Horne et al., 2009; Manz et al., 2020; Wang et al., 2019), including to predict mortality in patients with CKD (Forzley et al., 2018; Noh et al., 2020; Siga et al., 2020); it is not yet clear how such algorithms will compare to models such as ours that are based on clear theoretical constructs, or whether a hybrid approach may be most effective. Alternatively, EWSs based on clear theoretical constructs, both in medical practice and across disciplines, is largely based on level or dynamics of one or two variables that attempt to summarize the state of a complex system (Weinans et al., 2021). In the recent years, however, integrating larger numbers of variables into EWSs has started to gain attention (Dakos, 2018; Held and Kleinen, 2004; Lade and Gross, 2012; Lever et al., 2020; Suweis and D’Odorico, 2014; van de Leemput et al., 2014; Weinans et al., 2019). Our results demonstrate the potential of such multivariate integration, including with preexisting data in electronic clinical databases.

Here, our best models achieved an AUC of 0.82, a reasonable prediction given that all-cause mortality is a very broad and imperfect proxy for physiological collapse, and is notoriously challenging for clinicians to predict more than a few days in advance. In comparison, gait speed, a common test to screen for frailty and predict adverse outcomes in elders (Abellan Van Kan et al., 2009), generates AUC of 0.64 to 0.80 for predicting 5-year mortality in individuals aged 65 years or older, when combined to age and sex (Studenski et al., 2011). The increase in 3 months prior to death suggests a nearly ideal window for either clinical interventions to prevent decline, or for end-of-life discussions. However, it is not yet clear whether clinical interventions could successfully rescue patients from physiological collapse, though apparent spikes in CVPC1 before non-fatal hospitalizations (Figure 6) suggest hope in at least some cases. Alternatively, the best clinical application of EWSs may be in triggering discussions on end-of-life care preferences with patients. Access to quality palliative care and medical assistance in dying is becoming a priority (Government of Canada, 2020, 2018), and identifying the right time and way to initiate such a discussion still appears to be challenging for many nephrologists (Borreani and Miccinesi, 2008; Mandel et al., 2017). Once a sufficient algorithm and its application are validated, it should be possible to integrate into electronic health data systems (Major and Aphinyanaphongs, 2020).

Limitations of the study

While this study shows the clinical potential of biomarker variability as an EWS, more work is required before clinical implementation. Other multivariate indices of variability exist beyond CVPC1 (Liu et al., 2021; Weinans et al., 2021), and variability is only one signal among many that can be extracted from time series (e.g. lag-1 autocorrelations, flickering, skewness, etc., Scheffer et al., 2012). Future work will need to compare and likely integrate multiple approaches. Also, we could only include albumin in the 4-month list due to its measurement frequency in our study population, but our results suggest that if it were measured every two weeks it might improve the performance of integrative measures. Another important limitation of this study is the temporal resolution of 2 weeks; however, multivariate indicators based on variance perform better than those based on autocorrelation to predict critical transition with poor data resolution (Weinans et al., 2021). We included all-cause mortality since we did not have specific death causes and acknowledge that this is a limitation; nevertheless, our signal would likely have been stronger had we been able to specifically select mortality cases related to dialysis complications, given that unrelated mortality should be random and weaken the signal. Finally, our index appears to perform significantly less well when including the first six months on dialysis and patients who died within this period (HR95 of 3.2 compared to 9.7). These findings indicate that CVPC1 may not be the best indicator for short-term mortality after hemodialysis initiation, and might thus be preferred for long-term follow. Further studies should address this specific question, among others, before clinical implementation can be considered.

Conclusions

In conclusion, our findings have important implications at three levels. At a theoretical level, we showed surprising synchronicity in the variability of all measured biomarkers, implying that compartmentalization of physiology breaks down drastically preceding adverse critical transitions. At a methodological level, this suggests that the principles we illustrated can be leveraged to develop similarly powerful predictive indices in many other medical contexts—intensive care, diabetes, heart failure, and cancer, to name a few. At a practical clinical level, our findings imply that powerful predictors of mortality risk in patients with CKD under hemodialysis can be extracted from electronic medical record systems easily and within an appropriate time frame to help clinicians consider interventions or initiate a timely discussion with patients on end-of-life care.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information should be directed to and will be fulfilled by the lead contact, Alan A. Cohen (alan.cohen@USherbrooke.ca).

Materials availability

The study did not generate new unique reagents.

Experimental model and subject details

Our study population consisted of all patients treated at the hemodialysis external clinic of the CHUS hospital in Sherbrooke, Quebec, Canada between 1997 and 2017. CHUS protocol is hemodialysis 3x/week with a fixed blood work panel every two weeks to ensure follow-up and adjust treatment (Nakazato et al., 2017). Patient records were extracted from the CIRESSS (Centre informatisé de recherche évaluative en services et soins de santé) database, which aggregates all electronic hospital data for clinical and administrative purposes since 1991, a nearly exhaustive sample for the Eastern Townships region of Québec, population ∼326,000. This project received ethical approval by the Comité d'éthique de la recherche du CIUSSS de l’Estrie – CHUS (#2015-788, 14–059). Written informed consent for participation was not required for this study in accordance with the National Legislation and the Institutional Requirements, due to the retrospective nature of the study.

Method details

Cohort selection

To capture patients on whom sufficient data were gathered and with CKD rather than acute kidney failure, from 2,565 initial patients, we excluded 1694 patients who were no longer treated by in-center hemodialysis after six months (indicating death, recovery within six months from acute kidney failure or transition to another renal replacement modality or to conservative care), and 73 patients with irregular hemodialysis visits and an acute or unspecified kidney failure diagnosis. This left 798 patients with CKD on long-term hemodialysis. We excluded 28,681 visits (30%) with missing data (Figure 1 and Table S2). Because the start of hemodialysis itself may represent a critical transition (Broers et al., 2015; see also Figure S11), we also excluded data for the first six months on hemodialysis for each patient; however, for 26 patients, we were left with less than three visits and thus excluded these patients. Our final study cohort was comprised of 763 patients, for whom we included all visits for which biomarker data were gathered between earliest available hemodialysis visit (after the initial six-month exclusion) and death or last hemodialysis visit (Table 1). Our final study cohort included 479 (62.8%) men and 284 women. To validate our decision of excluding the first six months on hemodialysis, we also ran sensitivity analyses without excluding any observations, except the ones with missing data (n = 1765, Figures S5 and S6C), or by excluding only the first three months (n = 933, Figure S6B). Analyses were replicated in subpopulations to assess generalisability of findings: not diabetics (n = 372), diabetics (n = 391), men (n = 479), women (n = 284), <60 years old (n = 241), 60 to 75 years old (n = 368), and 75 + years old (n = 280). Age was defined based on age at visit, so some individuals were found in different categories over time. To assess if lost in follow-up impacted our findings, we ran sensitivity analyses excluding 112 patients who died but had no hemodialysis visit in the 30 days preceding their death (n = 651).

Missing data

To assess if data in our study cohort were missing completely at random, we performed Kruskal-Wallis test for age and Pearson's Chi-squared test for discrete covariates (sex, diabetes, and mortality), comparing observations with complete vs incomplete biomarker data. The significant differences in covariates suggest that our data are not missing completely at random (Table S2). Nevertheless, imputation was considered, but not used, for the following reasons:

Imputation methods have not been validated in this context, in contrast to more standard applications such as linear regression.

A successful imputation would have to preserve not only expected mean value, but expected variance of the ensemble, across relevant time periods. While it might be possible to develop methods that could achieve this, it would be a major undertaking in itself.

The creation of a potentially valid imputation method would require us to make numerous assumptions that might or might not be respected; as such, the risk of inducing a bias in the imputation is substantially larger than the risk of a bias due to the missing data.

If the result with imputed data conflicted with the current result, we would have no way of knowing for certain which method was preferable, but for reason 3) we would revert to the original analysis anyway.

The validation of the method and the stability of the results across numerous, distinct population subsets gives a strong indication that the result is not a function of population composition, which is the primary concern that would motivate imputation in this case.

Biomarkers and other key variables

We included all biomarkers measured regularly in the context of hemodialysis (Table 1), and excluded those measured only irregularly or having an important proportion of missingness in our cohort (uric acid, urea, iron binding capacity, iron, ferritin, iron saturation, transferrin, carbon dioxide, partial carbon dioxide pressure, partial oxygen pressure, pH, ionized calcium, intact parathyroid hormone, glycated hemoglobin and thyroid-stimulating hormone). For truncated laboratory variables (i.e. fields containing “<” or “>” signs), we generated the remaining part of the distribution (the tails of the distribution) using the rtnorm function from the MCMCglmm package (version 2.29, Hadfield, 2010). Mortality was identified via hospital and provincial (Régie d’assurances maladie du Québec) death records. Hospitalisation dates were recorded. Diabetes status was identified based on diagnostic codes 250 for ICD nine before April 2006, and E10, E11, E13, E14 for ICD-10 after April 2006.

Quantification and statistical analysis

Calculation of early warning signs

Biomarkers were log-transformed for normality as appropriate (Table 1). To test differences between biomarker levels and variability, and between univariate versus multivariate indices, we constructed four types of EWSs: two types of univariate indices and two types of multivariate (integrative) indices. Univariate indices were calculated either as within-individual means or CVs per time period, respectively measuring biomarker levels and variability (see Figure S13). Unless otherwise specified, CVs or means were calculated by year working backwards from time of death or censoring. We then combined information from individual biomarkers into multivariate indices using PCA applied on biomarker means or CVs. To calculate PCs based on biomarker means, we ran a PCA on all centered and scaled biomarker raw values using prcomp function in R, extracted the loadings, and applied them to the within-individual biomarker means calculated by year, thereby generating scores for each PC for every patient across years. CVs were calculated for each biomarker in each individual at each time period (by year, except for change point analyses) and were log-transformed. In a few cases, especially for tightly regulated biomarkers such as sodium and potassium, the CV per year was equal to zero because of laboratory rounding; for ease of analysis, we replaced these zeros by 90% of the next smallest CV value. We ran a PCA on the adjusted CVs (see below) using prcomp function in R, centering and scaling them. We noticed that CVs were biased by the number of values included in their calculation (Figure S13A). To confirm that this is a property of CV calculation, rather than of our data, we explored this association in simulated data. Simulated data were generated for each biomarker (transformed if needed to meet the assumption of normality) with the rnorm function, using the mean and standard deviation from our study population. Here, we only present results from one biomarker, but they all behave the same way: the fewer simulated values included in CV calculation, the lower CVs tended to be (Figure S13B). To control for this bias in CV calculation, we used the residuals from nonlinear regression models (CV=((1)⁄(√(n))×a+b), performed separately for each variable with the nls function in R). These residuals were used in the PCA to calculate CVPCs. We also ran sensitivity analyses using only CVs calculated with at least five observations, and using them directly in the PCA, instead of correcting with the non-linear regression (Figures S13C and S13D). These sensitivity analyses show that our results are strongly robust to these specific analytic details. Note that the square root of observation number included as weights in Cox proportional hazards models corrects for the loss of precision, whereas the non-linear regression described above corrects for the bias in CV calculation leading to smaller CV values when fewer observations are available, not for the precision.

Mortality prediction

Power of various EWSs to predict death was assessed using counting process Cox proportional hazards models (Andersen and Gill, 1982; Therneau and Grambsch, 2000). We used the hazard ratios to calculate “HR95” as a way to compare effect sizes across variables of varying scales: it represents the hazard ratio for an individual in the 97.fifth percentile of the variable relative to one in the 2.fifth, and permits apples-to-apples comparison across alternative biomarker indices on different scales (Milot et al., 2014). We assessed mortality risk with the coxph function (survival package version 3.1-12, Therneau, 2020), clustering multiple observations per individual using the “cluster” argument. Demographic covariates include age (modeled as a cubic spline with 5 degrees of freedom (df) using the bs function from the splines package to better account for non-linear relationships with age), sex, diabetes diagnosis, and years of follow-up. All models included the square root of observation number included in CVPC calculation as weights to account for the lower precision in CV estimation with fewer observations included in its calculation (note that the non-linear regression aforementioned does not correct for the loss of precision; it corrects for the bias in CV calculation leading to smaller CV values when fewer observations are available, not for the precision). We calculated the AUC with the pROC package (version 1.16.2) using a non-parametric method (Robin et al., 2011), while the Cox proportional hazards assumption was verified using the cox.zph function.

Change point analysis

We assessed the timing of changes in our integrative measures of biomarker variability (CVPCs) prior to death using change point regression with the mcp package (version 0.3.0, Lindeløv, 2020). Regression models were performed on all available CVPC values, and slopes were allowed to vary across individuals, with results presented after averaging values by time point (Figure 5).

Additional resources

Requests to access data: https://www.crchus.ca/en/services-outils/autresservices-et-outils/infocentre/.

Code: github.com/cohenaginglab/HD_variability.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Software and algorithms
R statistical language, version 3.6.0	The Comprehensive R Archive Network	https://www.r-project.org/
MCMCglmm package version 2.29	Hadfield, 2010	https://cran.r-project.org/web/packages/MCMCglmm/index.html
survival package version 3.1-12	Therneau, 2020	https://cran.r-project.org/web/packages/survival/index.html
pROC package (version 1.16.2)	Robin et al., 2011	https://cran.r-project.org/web/packages/pROC/index.html
mcp package (version 0.3.0)	Lindeløv, 2020	https://cran.r-project.org/web/packages/mcp/index.html