Identifying critical transitions in complex diseases.

Mortality and the burden of diseases worldwide continue to reach substantial numbers with societal development and urbanization. In the face of decline in human health, early detection of complex diseases is indispensable, albeit challenging. In this review, we document the research carried out thus far on the appearance of complex diseases marked by a critical transition or a sudden shift from a healthy state to a disease state. The theory of resilience and critical slowing down can provide practical tools to forecast the onset of various fatal and perpetuating diseases. However, critical transitions in diseases across diverse temporal and spatial scales may not always be preceded by critical slowing down. In this backdrop, an in-depth study of the underlying molecular mechanisms provides dynamic network biomarkers that can forecast potential critical transitions. We have put together the theory of complex diseases and resilience, and have discussed the need for advanced research in developing early warning signals in the field of medicine and health care. We conclude the review with a few open questions and prospects for research in this emerging field.

Introduction

The burden of diseases and associated fatalities on the global human population continues to rise due to environmental or behavioral attributes. Sudden deterioration in the health of individuals with the progression of diseases is common with hypertension (Lagro et al. 2012), diabetes (Li et al. 2014; Zeng et al. 2014), cancer of various organs (Kianercy et al. 2014), asthma attacks (Venegas et al. 2005; Winkler et al. 2015), epileptic seizures (Mormann et al. 2007), cardiovascular ailments (Kowalski et al. 2013), depression (van de Leemput et al. 2014), etc. Nevertheless, disease progression in a number of instances can be theoretically defined as a series of stages, viz., a healthy state (a stable state where the system encounters a gradual change on being perturbed), a pre-disease state (unstable or the edge of the healthy state before a drastic shift), and a stable disease state (Scheffer et al. 2009; Chen et al. 2012; Liu et al. 2012; Liu et al. 2013a, b). As the driver of disease increases beyond a critical threshold, the underlying system experiences an abrupt shift from a normal healthy state to a disease state. The risk accompanying the disease state of various complex diseases demands their early diagnosis, prevention, and cure. Detecting the pre-disease state is crucial to evade a critical transition to the disease state.

In the vicinity of a tipping point or a bifurcation point, critical transitions occur in the state of a system when tiny changes in an input condition cause a sudden shift to a contrasting state. For systems exhibiting fold or catastrophic bifurcation, a sudden shift to an alternate state is inevitable as the critical threshold for the control parameter (i.e., input condition) is crossed (Scheffer et al. 2001; Scheffer and Carpenter 2003; Scheffer 2009). A characteristic feature associated with such catastrophic events is the presence of the hysteresis loop (Scheffer et al. 2012). As the system collapses at the fold point with degrading conditions and in further trying to restore conditions by reversing the control parameter range, the point of recovery occurs past the point of collapse at another fold point. This phenomenon leads to the formation of a hysteresis loop. The distance between the fold points indicates the degree of hysteresis (Scheffer et al. 2001). Small changes in the vicinity of a tipping point can produce significant unforeseen changes in the state of a system (shift to a contrasting state). Anticipating such unannounced shifts in various real-world systems is challenging yet fundamental for the survival of the human race.

A plethora of mathematical methods incorporating stochasticity imitate the dynamics of complex systems exhibiting catastrophic bifurcation and multistability. Primitive theories confirm that a phenomenon commonly occurs in complex systems as they approach a tipping point – i.e., critical slowing down (CSD) (Tredicce et al. 2004). CSD is the phenomenon whereby the system’s rate of return to the current equilibrium state upon a random perturbation becomes increasingly slow as the dominant eigenvalue associated with the equilibrium state approaches zero. The autocorrelation at lag-1 reaches unity (Ives 1995; Dakos et al. 2008) and variance grows indefinitely (Carpenter and Brock 2006) as the eigenvalue goes to zero. Thus, an increasing trend in these statistical indicators works as an early warning signal (EWS) of an impending critical transition (Brock and Carpenter 2010; Dakos et al. 2012). Thereafter, researchers have developed generic as well as model-based EWSs and used them to forewarn of sudden shifts in diverse complex systems (Scheffer and Carpenter 2003; Lenton 2011; Veraart et al. 2012; Guttal et al. 2016). Recent studies have focused on transitions between alternate stable states in various biological systems (Lagro et al. 2014; van de Leemput et al. 2014; Sarkar et al. 2019). Slowing down has been used as a marker indicating the emergence of several chronic diseases (Trefois et al. 2015). Despite the usefulness of CSD-based generic indicators in forecasting critical transitions, they have their limitations (Ditlevsen and Johnsen 2010; Boettiger and Hastings 2012a, b; Clements et al. 2015; Dutta et al. 2018). This has led to the development of improved indicators such as dynamic network biomarkers (DNB) in preference over the traditional biomarkers (Chen et al. 2012; Liu et al. 2013a).

Recent studies have also detected critical transitions in various diseases using spatial patterns. Spatial patterns contain additional information depicting the system’s dynamics, providing predictions at a shorter lead time before the actual regime shift (van de Koppel et al. 2008; van de Leemput et al. 2015). Although less explored in disease biology, spatial early warning indicators have been used to forewarn of tippings in gene regulation systems (Yang et al. 2021), human seizures (Kramer et al. 2012), and various multiplex disease behaviors (Jentsch et al. 2018). The recent growth in the number of investigations on abrupt shifts in various diseases speaks of the importance of early diagnosis of the pre-disease state in the fields of medicine and health care. Albeit a well-developed theory exists on tipping points and their precursors, applying the latter in anticipating the fate of complex diseases is complicated. This entails a thorough understanding of the intertwined relationship among the various parameters triggering the system to a disease state.

In this review, we discuss sudden transitions observed in diverse systems of disease biology while deconstructing the vast literature on diagnosis and panacea for respective diseases. We study critical transitions in a two-dimensional cell-signaling system, the Cdc2-Cyclin B/Wee1 system which is a mutually inhibitory feedback loop (Sha et al. 2003; Angeli et al. 2004), in the presence of intrinsic noise. The activity and amount of the Cdc2-Cyclin B complex are associated with cells entering mitosis and causing cell proliferation. While unifying the mathematics associated with the critical transition in complex diseases, we also review the class of existing indicators. After that, we discuss the prospects, limitations, and future directions for this emerging area of research at the interface of complex systems, critical transitions, and disease biology.

Fundamentals of critical transitions

Resilience

In an ever-changing world, stochasticity is omnipresent. In multistable systems, due to the existence of multiple basins of attraction, small stochastic perturbations can allow a system to shift in the basin of another alternate state (Arnold et al. 1999). The probability or likelihood of the same is a function of the strength of perturbation as well as the size of the basin. Resilience is associated with the size of the valley/basin of attraction. It is a measure of the system’s ability to ingest the changes as the driver of decline increases, while ascertaining persistence; the deeper the valley, the larger is the perturbation it can withstand (figure 1c–d). The resilience of a system to random perturbations can be estimated by calculating extinction probabilities. This can be done by finding the domain of attraction regions where chances of shifts are more likely, displacing the system to the contrasting stable state. The approximate distance of the lowest point of the basin of attraction from the equilibrium point is a metric estimating the force necessary to drive the system to an extinction state (Holling 1973). In other words, resilience of a steady state depends on the distance to the boundary of the unstable manifold and how far away it is from the threshold value of the driver parameter (Scheffer et al. 2015; Strogatz 2018). The rate of recovery of the system (Pimm 1984) on being perturbed is estimated from the slope of the basin of attraction. Although the dominant eigenvalue approximates the slope, it does not completely determine the resilience of the system. As the system approaches the tipping point, CSD becomes more pronounced, indicating contraction of the basin of attraction of the present state (figure 1). This is accompanied by a decrease in recovery rate and decline in resilience. An alternative way to quantify the resiliance of a flickering system is by estimating the mean first-passage time (MFPT). MFPT from a basin of attraction can be calculated from the Fokker–Planck equation. However, estimating MFPT for real data is challenging as exit from a basin of attraction is hardly witnessed, apart from the need for massive data.

Figure 1

Critical transition in the concentration of protein Cdc2-cyclin B: (a) Stochastic time series of the model represented in equation 1 a-b: with feedback strength (far from the tipping point), and (b) with (close to the tipping point). (c) Schematic potential landscapes representing stable states of the deterministic system: higher resilience of the Cdc2-cyclin B state when it is far from the tipping point, and (d) lower resilience close to a tipping point, when the system approaches a sudden shift from the lower stable state to the upper stable state. (e, f) Closer to a tipping point, on account of reduced resilience, the system has more memory for perturbation than when far from a tipping point, characterized by higher SD and AR-1. (e) All other parameters for the circuit are given in the mathematical model section.

Critical slowing down

As stated in the introduction, CSD is a phenomenon where the system’s return to the current equilibrium state is slowed down upon perturbations in the vicinity of a tipping point. CSD occurs when the real part of the dominant eigenvalue approaches zero at the bifurcation point (Strogatz 2018). In general, the existence of a smooth bifurcation point in a dynamical system alone can determine the prevalence of CSD. The CSD property is prevalent in both catastrophic (fold) bifurcation and non-catastrophic (transcitical, pitchfork, or Hopf) bifurcation (Tredicce et al. 2004; Dutta et al. 2018). The dominant eigenvalue typifies the rate of change around equilibrium and is expected to essentially capture information related to tipping. Even though CSD may be present in systems far from the bifurcation point (Van Nes and Scheffer 2007), the property becomes more pronounced close to the latter as the recovery rate continuously decreases to zero. Slowing down energizes the system as the basin of attraction shrinks before eventually ceasing to exist. The rate of recovery on perturbing the system and the width of the basin of attraction are interrelated and vary with the timescale of the system dynamics. CSD observed in time series delineating a sequence of events prior to the aforementioned bifurcation forms the basis for development of early warning indicators of an impending transition (Scheffer et al. 2015). Theoretical studies confirm that dynamical systems with alternate stable states when observed at a spatial scale coupled through diffusion or advection of its state variables also exhibit CSD. It has been observed that CSD alters the diffusibility of the connected units in the system. This can be interpreted as the propensity of the system to return to its local equilibrium weakens as a critical threshold is reached (Dai et al. 2013). Taking into consideration spatial interactions, signatures of CSD can function to anticipate imminent transitions in a class of spatio-temporal systems (Dakos et al. 2010). While CSD may be persistent in a one-dimensional system or a network of interacting variables in spatio-temporal scales, it is not a precursor of all critical transitions (Dakos et al. 2015). Predictions using CSD-based indicators may be considered as the broader spectrum to only propel deeper investigation into the system dynamics and the change in resilience subject to infinitesimal perturbation close to equilibrium.

Sudden transitions in biological systems

Predictability of the fate of acute diseases depends on the knowledge of the complex nonlinear dynamics of specific disease systems. Although patterns of progression of diseases are diverse, there exist analogies of various human diseases with complex systems involving self-propagating positive and negative feedbacks among its state variables (Pomerening 2008). Being able to capture the stikingly similar dynamics of certain diseases, researchers have been able to identify abrupt shifts in chronic diseases and provide biomarkers forecasting their onset.

The genesis of several fatal diseases can be attributed to potential changes or degradation in gene expressions (Chen et al. 2012; Glass 2015; Trefois et al. 2015). Gene dysfunction involving the mutation, degradation, and modification of encoded proteins is majorly responsible for casualty due to hepatitis-B, hepatocellular carcinoma, and cancer of other organs (Yang et al. 2018). Recent studies have identified recurrent somatic mutation networks and they have performed analyses of gene mutations confirming its association with disease progression (Liu et al. 2013c). Positive feedback, a motif present in many genes found thus far, aids in the regulation of gene expression via synthesized proteins which promote the production of more proteins (figure 2). The presence of an auto-regulatory feedback module leads cells to exist in different steady states (Hari et al. 2020). Stochastic fluctuation in gene expression may trigger an abrupt shift from one state to another alternate state (Sharma et al. 2016; Sharma and Dutta 2017). Some instances are the existence of bistable gene expression in the induction of the lac operon in E. coli. An unannounced transition from an unregulated (state 1) to a regulated (state 2) state of the lac operon occurs at a critical value of an inducer concentration. In B. subtilis, positive feedback leads to the existence of low- and high-comM proteins. As the basal rate of protein synthesis crosses a threshold value, an abrupt shift is observed in the bacterial population level (Karmakar and Bose 2007). However, alterations in single gene expressions are often not sufficient to elucidate the course of action of a disease. Instead, gene regulation networks have reliably unveiled the progression of hepatitis-B (HBV)-hepatocellular carcinoma (HCC) and hepatitis-C (HCV)-HCC diseases (Yang et al. 2018).

Figure 2

Examples of circuits containing positive feedback loops: (a) Initiation of S-phase in eukaryote cell cycle. (b) Differentiation in ovaries of mammals. (c) Feedback in the two-component Cdc2-Cyclin B and Wee1 system. + sign denotes positive feedback.

Cancer involves abnormal cell growth and proliferation to other body parts. Almost of cancers of various organs trace their roots to dysregulation of pathways such as MAPK and P13K/AKT owing to multiple changes in several proteins (Samal et al. 2019). Although the origins of most solid tumors lie in epithelial cells, inhibitory feedback loops among gene expressions control epithelial–mesenchymal transitions (EMT). For instance, while the ZEB transcription factor family promotes EMT, the miR-200 family represses EMT. EMT and its reverse process,  mesenchymal–epithelial transition (MET), are important in cancer metastasis. Mathematical modeling of the ZEB/miR-200 feedback loop reveals multistability and critical transitions in cancer cells (Sarkar et al. 2019). Pulmonary disease such as asthma attacks also exhibit bistable regimes. Constriction of  bronchioles leads to reduced lung ventilation. Further increase in the critical driver value may lead to a catastrophic transition entailing severe breathing and ventilation defects (Venegas et al. 2005; Winkler et al. 2015). Congenital heart disease accounts for around of birth defect-related deaths. Unlike the adult circulatory system composed of atria and ventricles, the embryonic cardiovascular system consists of 6 pairs of aortic arches, out of which only 3 exist post stage-36. Investigation of intermediate stage-21 shed light on the critical transition in the embryonic circulatory system while estimating the susceptibility period (Kowalski et al. 2013).

Another deadly disease that perpetuates through an individual's lifetime is diabetes. There are two variants of diabetes: type 1 diabetes mellitus (T1DM) and type 2 diabetes mellitus (T2DM). In people suffering from T1DM diabetes, the person’s body is unable to produce insulin (essential for controlling blood sugar level), while in T2DM, one fails to adjust the blood sugar level even though insulin is produced. T1DM is a genetic disease which has a prolonged pre-diabetic period. Investigating the molecular mechanisms of the pre-diabetic state and distinguishing it from the other states is essential in evading a critical transition in T1DM (Liu et al. 2013b; Zeng et al. 2014). In T2DM, there exist 5 stages, viz., the latent stage, transition stage, impaired glucose tolerance (IGT) stage, impaired fasting glucose stage, and overt stage. While the patient is in any of the first 4 stages, recovery to the normal state is feasible. However, at the last stage, the system is unable to recover. In such a scenario, it is crucial to identify the transitions at the molecular level for early recovery using customized medication (Petersen et al. 2012). Loss of -cells that weaken connectivity of pancreatic islets act as a driver of critical transition (Szendroedi et al. 2007; Keller et al. 2008; Camastra et al. 2011).

The course of action of epidemic inflammatory diseases also exhibit bistable regimes where cytokinins act as a driver of inflammation (Ashley et al. 2012; O’Regan and Drake 2013). The emergence and spread of infectious diseases such as tuberculosis, measles, swine origin influenza A (H1N1), COVID-19, etc., are major health hazards affecting large parts of the community. An increase in (secondary infections from pre-infected individuals) over a threshold value can trigger a critical transition in a susceptible community. Researchers have studied both susceptible–infectious–susceptible (SIS) and susceptible–infectious–recovered (SIR), representing the dynamics of a large class of infections, by using master equations from statistical physics (O’Regan and Drake 2013; O’Dea and Drake 2019; Drake et al. 2019). Another fatal disease pathway that undergoes critical transition is epileptic seizure. This occurs due to malfunction of the nervous system accompanied by mild symptoms at an early stage to complete loss of control at later stages. Alternate stable states in seizure dynamics are the ictal and post-ictal states, as studied by Kramer et al. (2012) and Mormann et al. (2007). Recent studies have also found groups of neurons exhibiting Hopf bifurcation in seizures. Diseases such as depression, weakness in elderly people, and bacterial populations in the human intestine exhibit bistable steady states. Thus far, a wide range of diseases have been found to exhibit alternate stable states and a panacea to these complex diseases demand anticipating the pre-disease states and developing biomarkers that efficiently characterize their occurrence a priori. The phenomenon of CSD in temporal and spatial systems has been instrumental in developing genetic indicators that forewarn about transitions in many instances but does not always suffice. The complex dynamics of certain diseases require thorough investigation of the system and system-specific biomarkers.

Early warning indicators of sudden transitions in acute diseases

Early warning indicators for temporal systems

As a system approaches a bifurcation point, the return rate to equilibrium goes to zero. This is often expressed in terms of the dominant eigenvalue of the Jacobian matrix. From the fluctuation dissipation theorem, on approximating the critical point, variance grows large and autocorrelation at lag-1 reaches unity. Thus, CSD may be characterized by an increasing trend in variance and autocorrelation at lag-1 (figure 1e–f) closer to the tipping point, thus forewarning of a critical transition ahead of time (Ditlevsen and Johnsen 2010). As explained earlier, a change in the basin stability is observed as the system shifts to another steady state. This asymmetry leads to a skewed distribution of the stable states (Guttal and Jayaprakash 2008). Commonly used generic early warning indicators (Dai et al. 2012) of CSD are recovery rate (Wissel 1984; Drake and Griffen 2010), standard deviation (SD) (Carpenter and Brock 2006; Carpenter et al. 2008), autocorrelation at lag-1 (AR-1) (O’Regan and Drake 2013), skewness (Carpenter et al. 2008; Guttal and Jayaprakash 2008), kurtosis (Biggs et al. 2009), etc., and their spatial counterparts, viz., spatial correlation (Dakos et al. 2010, 2011), spatial variance (Guttal and Jayaprakash 2009; Dakos et al. 2011), etc. Trends in indicators of slowing down, viz., variance and temporal correlations across moods in healthy and depressed persons, are suggestive of using them as EWSs in mood tippings (van de Leemput et al. 2014; Olde Rikkert et al. 2016). Intuitively, for persons close to tipping (before the jump to depressed state), mood recovery post stress slows down, thus increasing AR-1 and variance. Slowing down also works as a metric for cardiac fitness tests (Olde Rikkert et al. 2016). Slow recovery of heart rate post exercise is an indicator of reduced fitness and may be a precursor to myocardial ischemia (Shephard 1967). Changing variance and skewness have proved to be potential indicators of critical transition in acute asthma attacks, forewarning of sluggish lung activity on ozone exposure (Venegas et al. 2005; Hsieh et al. 2014). Transitions to both emergence and elimination of such diseases have been anticipated using generic EWSs. These tools have been practical to an extent in forewarning about the outcomes of mitigation of COVID-19 and its impact on populations across countries (Kaur et al. 2020). Summary statistics such as CSD-based and likelihood-based tools for anticipating the emergence of diseases are more reliable when tested on high-resolution data. Detecting stronger signals for systems with fast timescales may require even more advanced methods (Brett et al. 2017). Critical transition in epileptic seizures are anticipated by rising AR-1 and flickering between ictal and post-ictal states (Kramer et al. 2012; Meisel et al. 2012). Flickering (Guttal and Jayaprakash 2008; Scheffer 2009) or bimodal frequency distribution have also signaled bistability in Cheyne–Stokes respiration and acute granulocyte leukemia (Olde Rikkert et al. 2016). For hibernating animals, flickering of body temperature signaled the transition between hibernation and activity states in species with gradual change in temperature (Oro and Freixas 2021). Frailty in individual patients remains understated by means of static metrics such as hospitalization, mortality, etc. However, dynamic resilience indicators such as AR-1, variance, and cross-correlation calculated on self-rated time series of individuals indicate frailty tippings in old age, which is further concerned with the resilience to various complex diseases (Gijzel et al. 2017). Apart from the diseases discussed above, recovery time predicts recovery post tippings in arrhythmia, colitis, falls, breast cancer (Stratton et al. 1982), neck cancer, smoking, and tube ventilation (Olde Rikkert and Melis 2019).

Early warning indicators for spatial systems

Theoretical studies have manifested the signatures of CSD to be also present in systems when observed at the spatial scale. Although envisioning systems across spatial scales aid in better understanding of the underlying mechanisms driving regime shifts, the lack of high-resolution spatial data and associated complexities have averted researchers from widely exploring spatial systems in various fields. However, in recent times, the availability of satellite and remote-sensing high-resolution data has kindled research of spatial systems (Kefi et al. 2014). Yang et al. (2021) studied critical transitions in a spatially extended gene regulation system where a sudden transition occurs from high to low protein concentrations. TF-A monomer concentration experiences a sudden collapse via a fold bifurcation as the transcription rate reaches a threshold value over time. This is supported by rising trends in variance, AR-1, and kurtosis at spatial scales. Asymmetry in spatial skewness is also observed prior to the transition in various dynamical systems (Yang et al. 2021). The various interventions to control protein concentration in the gene transcription regulation system have implications in providing remedies to a large class of complex diseases. Recently, the onset of epileptic seizure has been studied at different spatial scales (Meisel et al. 2012). The association of Hopf bifurcation with seizures and applicability of scaling laws at various spatial scales of a cluster of neurons and a rise in phase-locking near seizure thresholds have been reported. This has set goals to predict the smaller events, apart from the larger ones, leading to subclinical seizures. Another approach studied the dynamics of focal seizures and found that termination of seizures from the ictal to the post-ictal stage occurs via a discontinuous transition. These results are based on brain electrical activity data at different spatial scales. Although the termination of seizure was clearly confirmed by reduced neuron spiking, early warning indicators showed no trends. A rationale for this might be the microscopic dynamical mechanisms while the critical transition occur at a macroscopic spatial scale (Kramer et al. 2012). Thus, anticipating an impending transition in real-world datasets is often much more involved and requires better understanding of the spatial system and similar timescales. Studies also reveal the application of spatial indicators to multi-dimensional systems such as systems comprising the interaction of disease and public opinion. Tested across network types, increase in spatial correlation forewarn of critical transition from full vaccination (no disease state) to no vaccination (disease state). However, empirical social networks exhibit ambiguity in the type of transition, thus professing the need for further investigation to consider spatial CSD-based indicators for forewarning transitions in real networks (Dakos and Bascompte 2014). Although CSD-based indicators precede tippings in a range of complex diseases, a catastrophic transition in such diseases may occur without prior symptoms and early warnings (Dakos et al. 2015). However, system-dependent studies may better portray the course of the underlying dynamics of the disease while investigating stochasticity in molecular interactions, thereby paving the way for further research in developing biomarkers.

Dynamic network biomarkers

Traditional biomarkers were genes, molecules, or a network of genes that distinguish between a disease and normal state but fail to detect a pre-disease state due to their static property. However, this limitation has been overcome by the development of dynamic network biomarkers (DNBs), which can discern a pre-disease state from the normal state in complex diseases (figure 3). DNBs have successfully detected regime shifts in various biological processes related to the initiation of complex diseases, such as cell proliferation and differentiation (Zeng et al. 2013). As the system approaches a regime shift, groups of genes and transcription factors or molecules functioning as DNB exhibit increased correlation and considerable increase in variance among their expressions. DNBs for different diseases may not share common network topologies, although the statistical features such as SD and the Pearson correlation coefficient (PCC) remain the same, which aids in the identification of genes or cellular networks to function as a DNB. Based on the above criteria, a composite index is formulated that functions as an EWS for detecting a pre-disease state. These have been tested on experimental datasets of multi-period gene expression of three different tissues, liver, muscle, and adipose, in a diabetic rat model and control rat model. DNB methods not only provide early warnings for state transitions but also find the common function among these tissues from real spatio-temporal expression data (Li et al. 2014). However, the DNB for each tissue turns out to be different, stipulating the study of system-specific analysis. Time series of gene expression data in the mice metastasis model HCCLM3-RFP proclaims CALML3 as a DNB. The role of CALML3 in initiating metastasis was studied using two models – gain of function and loss of function. In the first model, HCCLML3 cells have low CALML3 concentration and high metastasis capacity. In the second model, HCCLML3 cells have high CALML3 concentration and low metastasis capacity. CALML3 acts as a biomarker anticipating pulmonary metastasis in liver cells, with loss in CALML3 providing indication for clinical treatment (Yang et al. 2018).

Figure 3

Schematic representation of the functioning of DNB (Dynamical Network Biomarker): Disease progression in several diseases occur in three states, i.e., the normal state, the pre-disease state, and the disease state. The normal state and the disease state exhibit high resilience, while the pre-disease state shows lower resilience.

Composite indices based on DNB have also proved fruitful in anticipating critical transitions in experimental data of acute lung injury, hepatitis-B, and lymphoma (Chen et al. 2017). An approach by Chen et al. (2017) considered temporal differential networks where edges connecting two nodes store correlation while edges connecting one node store variance. The inconsistency score (I-score) is obtained by training the hidden Markov model on a sequence of differential networks for N−1 time steps and then measuring the likelihood of a time point to be the end point of the process, and in this way marking the pre-disease state. The I-score effectively predicted pre-disease states in three real genomic datasets: acute lung injury, corneal trauma, and heregulin-induced breast cancer (Chen et al. 2017). The I-score provides robust prediction but is limited by the need for training with ample high-throughput data. Another interesting framework for EWSs for disease states is the progressive module network model (Zeng et al. 2014). This differs from classical DNBs by considering network constraints and is thus much more biologically relevant. The key features of the module include tissue and time specificity, and the network module is constructed based on a Markov cluster algorithm and from gene expression data, taking note of the possible biological interactions. This model has successfully detected the pre-disease state for T1DM from three different datasets having their origins in different parts of the body – pancreatic lymph nodes, spleen, and peripheral blood cells – based on the signal indices.

Predicting critical transitions in the Cdc2-Cyclin B/Wee1 system

Mathematical model

To better understand critical transitions and the applicability of EWSs, we present a detailed study of a cell signaling system involving proteins Cdc2-Cyclin B and Wee1. The Cdc2-Cyclin B/Wee1 system (see figure 2c) regulates mitotic control pathways and plays a significant role in cell division. Sudden changes in these protein concentrations regulated by the strength of the feedback loop can affect DNA repair and disease progression. As the feedback strength decreases, Cdc2-Cyclin B protein concentration increases, and further past a threshold value of v (see figure 4a), the system experiences an abrupt rise in Cdc2-Cyclin B concentration. High expression of Cdc2-Cyclin B can be involved in the process of tumor development and progression, which leads to disordered mitosis and malignant cell proliferation (Hou et al. 2018).

Figure 4

Critical transitions in protein Cdc2-Cyclin B and associated generic EWSs: (a) Transitions from a lower state to an upper stable state. Solid (cyan) lines indicate stable steady states, and dashed (red) lines indicate unstable steady states of the deterministic model. The black trajectory indicates stochastic time series. (b) Pre-transition stochastic time series (segment as indicated by the yellow boxed region in (a)). (c) Residual time series after applying a Gaussian filter (the orange curve in (b) is the trend used for filtering). (d and e) Generic EWSs calculated from the filtered time series after using a rolling window of of the data length: (d) variance and (e) AR-1. (f–i) Filtering bandwidth and rolling window chosen based on sensitivity analysis. (f) and (h) Contour plots showing the trends of generic EWSs variance and AR-1 respectively for different rolling-window sizes and filtering bandwidths as measured by the Kendall’s  value. The triangles indicate the rolling-window size and bandwidth used to calculate the EWSs. (g and i) Frequency distributions of Kendall’s  values corresponding to variance and AR-1, respectively.

The model system with a few underlying assumptions is given by equation (1), where x and y represent the protein concentrations of Cdc2-Cyclin B and Wee1, respectively (Angeli et al. 2004). The model can be written as: where , , , and  are rate constants. v represents the feedback strength varying in the range 0–2.  and are the Michaelis constants, and denote Hill’s coefficients. To study the influence of intrinsic noise, we derive the master equation for the grand probability function P(x, y, t) using the processes in table 1, which takes the form as below:where . The bifurcation diagram in figure 4a depicts the steady state density of protein-active Cdc2-Cyclin B as a function of the feedback strength (v). Cdc2-Cyclin B and Wee1 occur complementarily (a higher state for Wee1 implies a lower state for the other, and vice versa). A hysteresis loop occurs in the system on decreasing v from 2 to 0 and further increasing v in the opposite direction. The width of the hysteresis loop in marked by the bistable regions. We generated the time series under the influence of intrinsic stochastic fluctuations using Monte Carlo simulations. All the associated processes are presented in table 1. We calculated generic EWSs for the stochastic pre-transition time series, which may have implications for the early diagnosis of several acute diseases. We calculated the stochastic potential and mean first-passage time for the two-dimensional system to further analyze the dynamics before tipping.

Box 1

Glossary

Bifurcation: Gradual changes in the driver parameter causing a ‘qualitative’ change in the system dynamics.
Tipping point: A threshold value at which a dynamical system abruptly shifts from one stable state to another alternative stable state, in response to small stochastic perturbations.
Hysteresis: Multiple stable states exist under the same environmental conditions guided by a positive feedback loop. Reverting to the previous state by varying the bifurcation parameter in the opposite direction is non-trivial, i.e., point of collapse and recovery are not the same. This is particularly more common in systems with a saddle-node (fold) bifurcation.
Stochastic simulation algorithm: A Monte Carlo algorithm to generate numerical solutions of a stochastic equation system for known reaction rates. This is an exact scheme also known as the Gillespie algorithm.
Variance: The expectation of the squared deviation of a distribution about its mean. Slower recovery rates around the equilibrium on perturbations or flickering between alternate stable states increase the variance. Variance is defined as σz2=1n∑i=1n(z_i-μ)^2, where μ is the mean of the distribution.
Autocorrelation at lag-1 (AR-1): It is a measure of the change in short-term memory of a system prior to tipping. The formula to calculate AR-1 is given by: ρ=(E(z_t)-μ)(E(z_t+1)-μ)σz2. The slow rate of return to the equilibrium point on approaching the bifurcation point can be measured by change in the correlation of a time series.
Kendall’s τ: It indicates the strength of the trend in the generic indicators such as AR-1 and SD (σ_z) along the time series. The mathematical formula to calculate Kendall’s rank correlation τ statistics is as follows: τ=(c-d)/(c+d), where c is the concordant pair within the dataset, and d is the discordant pair within the dataset.
Stochastic potential: The negative logarithm of the probability of trajectories visiting a steady state. The probability is calculated as the ratio of frequency count of visiting a steady state to the total count of visits for a large number of initial conditions.
Mean first-passage time: The average time taken by the system to leave the basin of attraction. This is useful in characterizing the resilience of a system having alternate stable states or exhibiting flickering between states.

Table 1

Different birth and death processes, change of state vectors, gain and loss probabilities, and their propensity function, associated with the Cdc2-Cyclin B/Wee1 model represented by equation 1a–b

Sl. No.	Elementary events	Before reaction	After reaction	Gain probability	Loss probability	Propensity function
1.	ϕ→x	x-1y	xy	V_1α_1P(x-1,y)	V_1α_1P(x,y)	V_1α_1
2.	x+4y→z z→4y	x+1y	xy	β_1(vy)^4(x+1)k_1V14+(vy)^4P(x+1,y)	β_1(vy)^4xk_1V14+(vy)^4P(x,y)	β_1(vy)^4xk_1V14+(vy)^4
3.	x→ϕ	x+1y	xy	α_1(x+1)P(x+1,y)	α_1xP(x,y)	α_1x
4.	ϕ→ y	xy-1	xy	α_2V_1P(x,y-1)	α_2V_1P(x,y)	α_2V_1
5.	4x+y→z z→4x	xy+1	xy	β_2x^4(y+1)k_2V14+x^4P(x,y+1)	β_2x^4yk_2V14+x^4P(x,y)	β_2x^4yk_2V14+x^4
6.	y→ϕ	xy+1	xy	α_2yP(x,y+1)	α_2yP(x,y)	α_2y

 is the system’s volume where all the reactions occur. The symbols () and () in the column represent birth and death processes of the respective chemical species. Here, P stands for the grand probability function. and z are empty state and dummy variables, respectively.

Early warning signals

We calculated generic EWSs (Dakos et al. 2012) for anticipating abrupt shifts in Cdc2-Cyclin B concentration. Further, we comment on the reliability of the two most robust indicators, variance and AR-1, in forewarning about critical transitions in the above system. We calculated variance and AR-1 for the pre-transition stochastic time series (see figure 4). To avoid non-stationarities in time series that affect the prediction of generic EWSs, we applied Gaussian detrending using appropriate filtering bandwidth and rolling window to remove high frequencies and obtain the residuals (figure 4b–c). Rolling window and filtering bandwidth were chosen based on sensitivity analysis (see figure 4f–h) for variance and AR-1, respectively. As observed in figure 4d–e, both the indicators, variance and AR-1, successfully forewarn of an upcoming critical transition in the protein Cdc2-Cyclin B. Figure 4g–i shows a high Kendall’s  estimate for both the indicators. However, to comment on the reliability of the above indicators, we performed a surrogate analysis. We generated 1000 surrogate datasets and compared trend estimates to those of the original time series. We fitted the best linear auto-regressive moving-average model to the residuals for generating surrogates. We estimated the trends for variance and AR-1 using Kendall’s . Figure 5a–b shows that for the combination of window size and bandwidth  (of the original time series), AR-1 trends are not significant (p-value, 0.18), while trends for variance are significant with a p-value of 0.001. Thus, it might be worth highlighting that trends for CSD-based indicators may be considered, but with a note of caution, indicating the need for a further in-depth study of the respective system.

Figure 5

The probability distribution of Kendall’s  test statistic on a set of 1000 surrogate time series: Histograms depict the distribution of the test statistic for the surrogate time series (a) variance and (b) autocorrelation function at lag-1. Dashed (red) lines represent and confidence intervals. Solid (green) lines indicate the limit beyond which the Kendall’s  of the surrogate data is higher than the statistic observed in the indicators of the original time series. As observed, trends for variance are significant, while AR-1 trends are not significant.

Stochastic potential

We examined the resilience of the system to random perturbations by calculating the stochastic potential and basin stability (Sarkar et al. 2019) for the Cdc2-Cyclin B/Wee1 system (see figure 6). Higher values of potential indicate the existence of a shallow well, while lower values correspond to a deeper well. In other words, lower values in the color bar imply that on perturbing the system, the system has a higher probability of returning to the steady state. As shown in figure 6a–b, the potential well for the lower stable state corresponding to Cdc2-Cyclin B is deeper, while that for the higher stable state is shallower in the monostable regions (figure 6h–i). However, the depth of the potential for Cdc2-Cyclin B is greater for the upper stable state than the lower stable state. Consequently, a large perturbation is required to shift Cdc2-Cyclin B to the lower stable state as a deeper well signifies higher resilience to perturbations. Our conclusions based on the stochastic potential for the protein Cdc2-Cyclin B are also supported by the results obtained from calculating the basin stability for the system. For a sufficiently large number of random initial conditions, basin stability is the probability of returning to the original state. In the monostable regions, this probability stands at 1 (see figure 6c–j). However, in the bistable regions, for different initial conditions, trajectories reach the upper stable state with probability 0.75, while they shift to the lower stable state with probability 0.25 (figure 6g). This can be attributed to the deeper well of the upper stable state of the protein Cdc2-Cyclin B in the bistable regions as indicated by the stochastic potentials in figure 6f.

Figure 6

Stochastic potential landscapes and basin stability for feedback strength (v) obtained from the master equation (2): Stochastic potential for: (a) monostable low-density state for Cdc2-Cyclin B at , (e) bistable high-density and low-density Cdc2-Cyclin B states for , and (h) monostable high-density Cdc2-Cyclin B state at . The blowup diagrams (a), (d), (f) and (i) represent magnified regions in potential wells. The color bar represents the negative logarithm of steady state probability [i.e., ]. (c, g and j) Pie diagrams representing basin stability of the system for feedback strengths , 0.9, 0.4, respectively. The basin stability measure is calculated for percentages of simulations with random initial conditions reaching a particular steady state in a monostable or bistable region. Green and peach regions correspond to the percentage of simulations reaching upper and lower states, respectively. x and y denote Cdc2-Cyclin B complex and Wee 1, respectively.

Mean first-passage time of the Cdc2-Cyclin B/Wee1 system

We explored the average time taken by the protein Cdc2-Cyclin B to reach an alternate stable state for the first time, having started at an initial stable state (Samal et al. 2021). For systems existing in the bistable region, the first-passage time (FPT) gives an estimate of the time taken by the system to leave a potential well for a chosen initial point for the first time. MFPT is the average of the FPTs. We derived MFPT for the Cdc2-Cyclin B state variable numerically using the Gillespie simulations for a set of initial conditions in a potential of volume V; we noted the time at which the state variable exits V for the first time. The inverse of MFPT quantifies the rate at which the state variable Cdc2-Cyclin B hits the potential barrier (unstable steady state). MFPT for Cdc2-Cyclin B from upper to lower states and the converse are reported as the average of 5000 realizations sampled all over the potential well, and all the initial conditions are taken uniformly from the ranges [0, 1.2] × [0, 1.2]. As observed in figure 7, for both the transitions, as we move closer to the bifurcation point, the time taken to shift to an alternate state for the first time decreases. Evidently, as the feedback strength decreases, the concentration of Cdc2-Cyclin B increases, and at the threshold value of v, the system jumps to an alternate state. This increased amount causes cells to enter mitosis, losing control over cell division and promoting cancer cell proliferation and development conditions. Hence, MFPT provides us insight into the time before which necessary control measures can prevent the abrupt shift of the Cdc2 complex.

Figure 7

MFPT with decrease in feedback strength (v): The red curve represents the MFPT calculated for Cdc2-Cyclin B system to switch from the upper state to the lower steady state under the influence of intrinsic noise. Similarly, the blue curve shows the MFPT to switch from the lower steady state to the upper state.

Prospects

CSD can form the basis for anticipating critical transitions associated with complex diseases a priori, but considering it a panacea for forewarning all sudden shifts is misleading. Moreover, the CSD theory is well established for mathematical models representing complex stochastic systems, but their applicability to real empirical data remains limited. For various dynamical systems in bistable regions perturbed by external noise, the efficacy of generic EWSs can vary as shown in the literature (Qin and Tang 2018; Samal et al. 2021). In fact, for external noise, the efficacy of EWSs can reduce, and the effect become more profound with increasing noise amplitude and correlation timescale (Qin and Tang 2018). In instances where CSD-based indicators provide early warnings, the signals should be relied upon for further considering the system for an in-depth study of the molecular mechanisms that guide it en route to the alternate stable state. Studying the dynamics of such models and dominant eigenvalues can aid in CSD prediction. However, to mitigate a critical transition in cancer patients or other diseases, it is important to identify the bifurcation parameter (small gradients of the parameter which can produce large changes on crossing a threshold value) (Korolev et al. 2014). An in-depth study of the evolution in human cancer samples can aid in more precise modeling of tumor dynamics. However, despite the complexities, researchers have put significant effort into developing measurable signals of tipping points in disease biology. This has led to the development of network biomarkers, including single-gene biomarkers, or involving groups of genes, and further advancement has brought DNBs into the field. Although DNBs for systems decipher information at the molecular level, seeking a reliable indicator out of them requires high-throughput empirical data. Nevertheless, the success of DNBs is due to the fact that changes in the system dynamics near a critical threshold may be envisioned through changes in gene (protein) sequences. At the same time, generic CSD-based indicators only capture fragments of it. A key goal in the health-care sector is to increase adaptability and recovery on exposure to stressors. Achieving this requires identifying drivers and their interconnections while developing improved DNBs.

Building resilience is an ongoing process and shall continue at various levels of health care with the emergence of complex diseases (Martin 2018). As conservationists, we need to learn from the prevailing conditions to forewarn and hence be fortified (Scheffer et al. 2018). A potential advancement in the field to achieve resilience is the use of machine learning (ML). ML models can be trained using data from an array of critical behaviors in biological systems to predict the onset of critical transition. While ML methods have shown their mystique in discerning transitions (Dev et al. 2021), it remains to apply ML techniques to disease transitions. Lastly, developing robust and reliable indicators to anticipate instabilities in a ceaselessly changing world requires real-time monitoring of systems rather than predicting in retrospect. This upfront area of research leaves open questions related to sudden transitions in disease biology. A few major questions that need to be addressed are: Is it possible to devise indicators that optimize the timescale while providing reliable early warnings before a disease state? Can we develop an absolute measure of resilience in biological systems amidst noisy sequences of intertwined events?