 Research
 Open Access
The decrease of consistence probability: at the crossroad of catastrophic transition of a biological system
 Pei Chen^{1} and
 Yongjun Li^{1}Email author
https://doi.org/10.1186/s129180160295y
© The Author(s) 2016
Published: 1 August 2016
Abstract
Background
Unlike traditional detection of a disease state in which there are clear phenomena, it is usually a challenge to identify the predisease state during the progression of a complex disease just before the serious deterioration, not only because of the high complexity of the biological system, but there may be few clues and apparent changes appearing until the catastrophic critical transition occurs.
Results
In this work, by exploiting the different dynamical features between the normal and predisease states, we present a hiddenMarkovmodel (HMM) based computational method to identify the predisease state and elucidate the essential mechanisms during the critical transition at the network level. Specifically, by considering the network variation and regarding that the predisease state is the end or shiftpoint of a stationary Markov process, a consistence score is proposed to measure the probability that a system is in consistency with the normal state. As validation, this approach is applied to detect the upcoming critical transition of complex systems based on both the dataset generated from a simulated network and the rich information provided by highthroughput microarray data. The effectiveness of our method has been demonstrated by the identification of the predisease states for two real datasets including HCVinduced hepatocellular carcinoma and virusinduced influenza infection.
Conclusion
From dynamical view point, the criticaltransition phenomena in many biological processes are of some generic properties, which can be detected by the established method.
Keywords
 Dynamical network biomarker
 Hidden Markov process
 Predisease states
Background
Recently, evidence suggests that the deterioration of many complex diseases is not necessarily smooth but abrupt, that is, the sudden change of system state exists widely during the progression of complex diseases. For example, some chronic diseases such as cancer, the malignant deterioration may arise within a period of shorttime progression, while before such catastrophic transitions the disease such as chronic inflammation may progress gradually for years of long incubative duration [1–5]. In other words, during the progression of illness there is a sudden critical state transition from a relatively healthy stage to a seriously diseased stage. For many complex diseases, it is crucial to detect such critical state transition in advance so as to prevent or at least get ready for such a catastrophic event. However, it is still a challenge work to signal the upcoming critical deterioration since the state of the system may show little apparent change before the tipping point is really reached. This is also the reason why diagnosis based on traditional biomarkers may fail to indicate a predisease state. A possible approach to study the warning signal of the sudden deterioration is to explore and analyze the dynamical features generated from the early abnormalities in distinct timeseries prior to the emergence of the apparent malignancy. Therefore, in order to describe the underlying dynamical mechanism of complex diseases, their evolutions are often modeled as timedependent nonlinear dynamical systems, in which the abrupt deterioration or qualitative transition is viewed as the state transition or phase shift at a bifurcation point [6]. We particularly focus on the complex diseases with sudden deterioration phases or critical transition points during their progressions.
In this work, by exploring the distinct dynamical features between the correlation networks respectively generated in normal and predisease state, we developed a computational method on the basis of the hidden Markov model (HMM) for identifying the predisease state before the critical point is really reached during the biological process of complex diseases. Specifically, it is natural to model the progression of a biological system in a normal state as a stationary Markov process, since the normal state is a stable state and with high resilience. The predisease state is modelled as the timevarying Markov process due to its unstable nature and high sensitivity to even small perturbation. The disease state is another stationary Markov process in view of its high stability (see Fig. 1 a). Identifying the predisease state is then equivalent to detecting the end of the stationary Markov process. Utilizing the timecourse data, we presented the computational method and algorithm on estimating the possibility of supposed termination of Markov process at each candidate sampling point. Specifically, by exploring the critical phenomena of network structure in dynamics (Fig. 1 b), a consistence score (Cscore) was proposed to signal the upcoming critical transition, i.e., the drastic decrease of Cscore implies the onset of a predisease state, in contrast to the relatively smooth Cscore in either a normal or disease state (Fig. 1 c). To demonstrate the effectiveness of our method, we applied the algorithm to a simulated regulation network and two sets of real data, the microarray dataset of HCVinduced dysplasia and hepatocellular carcinoma (HCC) (GSE6764) and live influenza infection (humans) caused by H3N2 virus (GSE30550). The predisease states were successfully identified for both numerical simulation and real datasets, and thus signaling the imminent critical transitions.
Methods
We first present the theoretical basis, i.e., the dynamical properties of a complex system near the tipping point, and then illustrate the preprocessing of real datasets and the detail algorithm.
Theoretical basis
Disease progression or its biological process can be generally divided into three states or stages, i.e., (A) the normal stage, (B) the predisease stage, and (C) the disease stage (Fig. 1 a). The normal stage is a stable state with high resilience and robustness stage, during which the state may change slowly and thus is modelled as a stationary Markov process. The predisease stage is unstable and defined as the limit of the normal stage just before the occurrence of catastrophic phase shift. It is sensitive to perturbation including noise or external interference that leads to the change of system parameters, thus still reversible to the normal stage given appropriate interventions. Therefore, the system progression during a predisease stage is considered as a timevarying Markov process, during which the statetransition probability may fluctuate from time to time. However, further progression of the illness led by persistent effects of perturbation may trigger a drastic state change into the disease stage, the other stable state described as the second stationary Markov process, which is usually difficult to return to the normal state even with intensive interventions. Hence, it is crucial to detect the predisease state so as to prevent qualitative deterioration into an irreversible stage. On the basis of the above settings, detecting the imminent critical transition is equivalent to identifying the end point (or switching point) of the stationary Markov process (Fig. 1). Besides, we investigate the different dynamical features between the correlation network respectively generated from normal and predisease state, i.e., comparing the differential links from adjacent time points.
Based on such study design, we carry out theoretical derivation in the following sections.
Markov process of the network evolution near the critical point
We describe the theoretical derivation of our computational method, and introduce the qualitative behaviors in dynamics of biological variables to characterize the critical transition. The dynamics for the progression of complex diseases is very complicated either before or after the critical transition, and therefore the state equations are generally constructed in a highdimensional space with a large number of variables and parameters. Therefore, it is a difficult task to construct an accurate mathematical model describing the dynamical behavior of the system during the biological process. Thus we aim at developing a modelfree method to detect the critical signal.
where Z(k)=(z _{1}(k),…,z _{ n }(k)) is an ndimensional state vector or variables at time instant k that represents gene or protein expressions, while P=(p _{1},…,p _{ s }) is a parameter vector or driving factors that represent slowly changing factors, e.g., genetic factors (SNP, CNV, etc.) and epigenetic factors (methylation, acetylation, etc.). f:R ^{ n }×R ^{ s }→R ^{ n } are generally nonlinear functions. Furthermore, the following conditions are assumed to be held for system (1). (1) \(\bar {Z}\) is a fixed point of system (1) such that \(\bar {Z}= f(\bar {Z}; P)\). (2) There is a value P _{ c } such that one or a pair of eigenvalues of the Jacobian matrix \(\left.\frac {\partial f(Z; P_{c})}{\partial Z}\right _{Z=\bar {Z}}\) is equal to 1 in the modulus. (3) When P≠P _{ c }, the eigenvalues of (1) are not always equal to 1 in the modulus. These three conditions with other transversal conditions imply that the system undergoes a phase change at \(\bar {Z}\) or a codimensionone bifurcation when P reaches the threshold P _{ c }.
For system (1) near \(\bar {Z}\), before P reaches P _{ c }, the system is supposed to stay at a stable fixed point \(\bar {Z}\) and therefore all the eigenvalues are within (0,1) in modulus. The parameter value P _{ c } at which the state shift of the system occurs is called a bifurcation parameter value, or a critical transition value.
where ζ=(ζ _{1},…,ζ _{ n }) are small Gaussian noise with zero means. ζ _{ i } has a small standard deviation σ _{ i } for all i, and covariances κ _{ ij }=Cov(ζ _{ i },ζ _{ j }).

when s _{ i1}≠0 and s _{ j1}≠0, \(\lim \limits _{\lambda _{1}\rightarrow 1}\text {PCC}(z_{i}, z_{j})\rightarrow 1\);

when s _{ i1}≠0 and s _{ j1}=0, \(\lim \limits _{\lambda _{1}\rightarrow 1}\text {PCC}(z_{i}, z_{j})\rightarrow 0\);

when s _{ i1}=0 and s _{ j1}=0, \(\lim \limits _{\lambda _{1}\rightarrow 1}\text {PCC}(z_{i}, z_{j})\rightarrow P_{ij}\), where P _{ ij } is a bounded value.
Hence, close to a tipping point, among the original variables Z=(z _{1},...,z _{ n }) there is a dominant group which is composed of dominant variables z _{ i }=s _{ i1} y _{1}(k)+⋯+s _{ in } y _{ n }(k) with s _{ i1}≠0. It is clear from the above derivation that the correlation between a pair of dominant variables increases sharply as the dominant eigenvalue λ _{1}→1, while the correlation between a dominant variable and any other molecule decreases sharply. It also should be noted that such critical change of the correlation only appears when the system approaches to the critical tipping point, or equivalently, the system is in a predisease state. By employing this dynamical feature between a normal state (when the system is far from the tipping point) and a predisease state (when the system is in the vicinity of the tipping point), it is possible to detect the earlywarning signal of the critical transition based on the hidden Markov model.
Identifying the end of Markov process and the algorithm of HMMbased method

Denote the stationary Markov process as M _{1}, and the timevarying Markov process as M _{2}.

Denote the time variable as t, and the progression of the system along time series as t∈{1,2,...,T−1,T,...}.

Denote the observed sequence up to time point t as O={o _{1},o _{2},...,o _{ t−1},o _{ t }}, where o _{ t } represents the sample set derived at time point t.

Denote the state sequence up to time point T as {s _{1},s _{2},...,s _{ T−1},s _{ T }}, i.e., the state of the system is s _{ T } at time point t=T, or equivalently, s _{ T }=S t a t e(o _{ T }).
For each candidate time point t=T, the high value of Cscore presents that the progression of the system at t=T is consistent with the stationary Markov process M _{1}, i.e., it is still in the stationary Markov process, while the sudden decrease of Cscore illustrates the low consistence with M _{1} (Fig. 1 c), and the progression of the system is no longer in the stationary Markov process. Therefore, the abrupt change of Cscore identifies the predisease state and indicates the upcoming critical transition.
 1.
Train a hidden Markov model (HMM) θ _{ T−1}=(A,B,π) on the basis of an observed sequence {o _{1},o _{2},...,o _{ T−1}}, i.e., the preceding T−1 sets of samples generated from time points 1,2,...,T−1. The stationary Markov process in the normal state is actually described by the trained HMM.
 2.
Calculating the Cscore based on the observation {o _{ T }} and the trained HMM θ _{ T−1}. If there is a drastic decrease of Cscore, then the iterative process end up with t=T being the switching point, at which the biological system is in the predisease stage. Otherwise go to back to the training step for next time point t=T+1.
First, to train an HMM θ _{ T−1}=(A,B,π) where the subscript T−1 of θ represents that the HMM is derived from the training samples up to time point t=T−1, we need to estimate a state transition matrix A, an emission matrix B, and a probability vector for the initial state π. For a network with n nodes and m links where each node represents a biomolecule and each link represents the correlation between two nodes, suppose at a sampling time point t∈{1,2,...,T−1,T,...} there are w samples for each node z _{ i }, i.e., \(\{{z_{i}^{1}} (T1),{z_{i}^{2}} (T1),...,{z_{i}^{w}} (T1)\}\). Then through leavingoneout procedure we obtain w Pearson’s correlation coefficients (PCCs) between any two nodes z _{ i } and z _{ j }, i.e., {P C C _{1}(z _{ i },z _{ j }),P C C _{2}(z _{ i },z _{ j }),...,P C C _{ w }(z _{ i },z _{ j })} where each P C C _{ k }(z _{ i },z _{ j }) (k=1,2,...,w) is calculated based on w−1 samples for z _{ i } and z _{ j }. To train A and B based on an unsupervised learning procedure, we have the following steps.
A. Estimate the distribution of each link at a former time point (T−2) Under the assumption that each correlation coefficient follows Gaussian distribution, we obtain the estimation of the distribution for link P C C(z _{ i },z _{ j })_{ T−2} between two nodes z _{ i } and z _{ j } at time point T−2, i.e., based on the w correlation coefficients, we estimate the mean μ _{ k }(T−2) and standard deviation σ _{ k }(T−2) for each link P C C(z _{ i },z _{ j })_{ T−2}. Then we have the distribution \(N(\mu _{k} (T2), {\sigma ^{2}_{k}} (T2))\).
Obviously, \({L_{k}^{s}} (T1)=0\) represents that the correlation l i n k _{ k }(T−1) is consistent with the former distribution \(N(\mu _{k} (T2), {\sigma ^{2}_{k}} (T2))\), while x _{ k }(T−1)=1 represents that the correlation l i n k _{ k }(T−1) is inconsistent with the former distribution \(N(\mu _{k} (T2), {\sigma ^{2}_{k}} (T2))\). Thus, for each sample of correlation \((lin{k_{1}^{s}} (T1),lin{k_{2}^{s}} (T1),...,lin{k_{m}^{s}} (T1))\), the vector \(L^{s} (t1)=({L_{1}^{s}} (T1),..., {L_{m}^{s}} (T1))\) is the consistence vector at time T−1.
Let #0(T−1) and #1(T−1) respectively denote the number of value 0 and that of value 1 in an consistence vector L ^{ s }(T−1) at T−1. Obviously, #0(T−1) + #1(T−1)=m, where m is the number of links in the network, among which there are #0(T−1) variables consistent with the former distribution \(N(\mu _{k} (T2), {\sigma ^{2}_{k}} (T2))\), while #1(T−1) variables inconsistent with the former distribution \(N(\mu _{k} (T2), {\sigma ^{2}_{k}} (T2))\).
According to above settings, we actually transform the observed correlation sample set o _{ T−1}=(l i n k _{1}(T−1),l i n k _{2}(T−1),...,l i n k _{ m }(T−1)) into the corresponding consistence vector o _{ T−1}=(L ^{1}(T−1),L ^{2}(T−1),...,L ^{ m }(T−1)).
C. Training the HMM at T−1 In this step, we need to identify the state transition matrix A and the emission matrix B at (T−1), that is, training the HMM θ _{ T−1}=(A(T−1),B(T−1),π) on the basis of an observed sequence {o _{1},o _{2},...,o _{ T−1}}.
with i,j∈{1,2}.
where j∈{1,2} and k∈{0,1,2,...,m}. Obviously, there are m+1 possible observable cases for any correlation sample at t−1, i.e., case #1(T−1)=k with k∈{0,1,2,...,m}. In the case of an mlink biological network, case #1(T−1)=k reflects that there are k links differentially expressed in one observation (i.e., one sample) at T−1 comparing with their former expressions.
with i∈{1,2}.

Initialization For h=0, set initial values for \(a_{ij}^{0}\), \(\,b_{jk}^{0}\), and \({\pi _{i}^{0}}\), we have the HMM θ ^{0}=(A ^{0},B ^{0},π ^{0}).

Update For h=1,2,..., we have the update for \(a_{ij}^{h}\), \(\,b_{jk}^{h}\), and \({\pi _{i}^{h}}\) by recursion$${} a_{ij}^{h} =\frac{\sum\limits_{t=1}^{T1}\xi_{t}(i,j)}{\sum\limits_{t=1}^{T1}\gamma_{t}(i)}, \quad \!\!\!b_{jk}^{h} =\frac{\sum\limits_{t=1, \#1(T1)=k}^{T1}\gamma_{t}(k)}{\sum\limits_{t=1}^{T1}\gamma_{t}(k)}, \quad{\pi_{i}^{h}} =\gamma_{1}(i), $$(9)where$$ \gamma_{t}(i)=P(s_{t}=M_{i}\,\,O,\, \theta_{p})=\frac{P\left(s_{t}=M_{i},\, O\,\,\theta_{p}\right)}{P(O\,\,\theta_{p})} $$(10)and$$ {} {{\begin{aligned} \xi_{t}(i,j)=P\left(s_{t1}=M_{i}, s_{t}=M_{j}\,\,O,\, \theta_{p}\right)=\frac{P\left(s_{t1}=M_{i}, s_{t}=M_{j},\, O\,\,\theta_{p}\right)}{P\left(O\,\,\theta_{p}\right)} \end{aligned}}} $$(11)
with i,j∈{0,1}. For γ _{ t }(i) and ξ _{ t }(i,j), the HMM θ _{ p } used in the prior knowledge is that updated from the preceding step. For example, at the first iterative step, the HMM θ _{ p } is θ ^{0}=(A ^{0},B ^{0},π ^{0}) based on the initial values. The observation sequence used in the prior knowledge is O={o _{1},o _{2},...,o _{ T−1}}.

Ending When h=H, i.e., the Hthupdating step, the recursion is terminated. Then$$ {\theta_{i}^{H}} =\left(A^{H},B^{H},\pi^{H}\right). $$(12)
The HMM used in the testing process follows \(\theta _{T1}={\theta _{i}^{H}}\).
Under the assumption that the transition point is at T, or in other word, time point T is hypothesized as the end point of a stationary Markov process of the normal stage (see Fig. 1). Thus the onset of a predisease stage is the end of the stationary Markov process described as the trained HMM θ _{ T−1}. Therefore, at testing step in a candidate transition point T, we calculate the consistence score, i.e., Cscore, based on the trained HMM θ _{ T−1}=(A(T−1),B(T−1),π).
where a _{11} and a _{ ij } is from the state transition matrix A=(a _{ ij })_{2×2} in Eq. (6), b _{1k } and b _{ jk } is from the emission matrix B=(b _{ jk })_{2×(m+1)} in Eq. (7) while k=#1(T) represents that for the sample set o _{ T } there are k variables with consistence index 1 in average, Q is the forward probability calculated based on standard forward algorithm. It should be noticed that in Eqs. (13) and (14) the backward probability is set to be 1, since samples o _{ T+1},⋯ are not available when T is the testing time point.
According to above settings, given the HMM θ _{ T−1}, the calculation of HMM probability P _{ T } at a candidate time point T only relies on the samples from T−1 and T. Obtaining the Cscore P _{ t } for every candidate time point, the time point a r g _{ t }[m a x(P _{ t })]_{ t=1,2,...,T }, is the transition point.
Results
Identifying the pretransition state for a sevennode network
where q is a scalar control parameter and ζ _{ i }(t) (i=1,2,...,10) are Gaussian noises with zero means and covariances κ _{ ij }=Cov(ζ _{ i },ζ _{ j }). z _{ i } (i=1,...,10) represent the concentrations of mRNAi. In Eq.(15), the degradation rates of mRNAs are \(\left (\frac {2+q}{5}, \frac {2+q}{5}, 1, \frac {6}{5}, \frac {7}{5}, \frac {9}{5}, \frac {9}{5}\right)\). There is a stable equilibrium point \(\bar {Z}=(\bar {z_{1}}, \bar {z_{2}},...,\bar {z_{10}})=(0,0, 0, 0, 0, 0, 0, 0, 0, 0)\). The differential equations Eq. (15) can be transformed into the difference equations Z(k+1)=f(Z(k),P) using Euler scheme with a small time interval 1. It is clear that there are seven distinct eigenvalues (0.67^{q},0.45,0.37,0.30,0.24,0.20,0.13) for the linearized system. Thus, the equilibrium point \(\bar {Z}\) is stable when q∈(0,1]. There is a critical value q _{ c }=0. We aim to detect early warning signals that indicate the critical transition as a control parameter q approaches the critical value 0 from q>0. Applying the HMMbased approach to the system, we obtain the Cscore curve as in Fig. 2 b.
The numerical simulation shows that a drastic boost of the Cscore, i.e., HMM probability, indicates the upcoming critical transition at parameter q=0 (Fig. 2 b). To demonstrate the different dynamics of the system between the normal state and the predisease state, we illustrate the underlying frequency of links with different correlation values (Fig. 2 c), from which it can be seen that there is a significant change in the frequency distribution of the links when the system is near a tipping point.
Predicting critical transitions in real datasets
We applied the HMMbased method in three real experimental datasets, i.e., the microarray data for HCVinduced dysplasia and hepatocellular carcinoma (HCC) (GSE6764) and live influenza infection (humans) caused by H3N2 virus (GSE30550).
The functional analysis shows that some of the selected genes are highly relevant to the corresponding complex diseases, which validates the effectiveness of our method in a way. In the HCC study, many genes included in the top significant subnetwork relate to the response to HCV infection, especially the activation of the immune system and the dysfunctions associated with basic cell metabolism of hosts [23–25]. In the enrichment analysis, the most significant enriched pathways are related to the function of cell growth and cell metabolism, such as transcriptional misregulation in cancer, the Wnt signaling pathway and purine metabolism. The enriched pathways in cancer and hepatitis provide evidence that most of the genes related to differentiallyexpressed links may relate to the deterioration into HCC. These functional analysis implies the involvement of the differentialexpression genes and links in the dysfunctional pathways and other HCVrelated biological processes.
Figures 3 b, 3 d and 4 b shows another application of Cscore in the dataset of H3N2 virusinduced influenza infection, in which there are 16 sampling time points over the whole study period (132 hours). Nine subjects were diagnosed as having influenza infection or corresponding clinic symptoms 45 hours after the exposure to influenza viruses [26]. The specific procedure of data processing, gene filtering and computation are similar to the previous application. It can be seen that the Cscore curves based on the human PPI network for live influenza infection in Fig. 3 b, with eight probability curves respectively corresponding to the first eight candidate points. Among the probability curves in Fig. 3 b, the red curve presents the Cscore based on the top 5 % differentialexpression genes at 29 hr, while the orange one shows that calculated at 36 hr, the adjacent time point after 29 hr. Both these two curves show a sudden decrease of consistence probability during the progression, which implies that the onset of predisease state is around a period spanned from 29 to 36 hr, i.e., the upcoming deterioration into a disease state might be after 36 hr, which is in coincidence with the fact that the early symptoms of influenza infection arises after 45 hr. Furthermore, to show the significance of the selected genes whose collective dynamics results in the significant changes in the links and thus generate the earliest signal at 29 hr, in Fig. 3 d we carried out the clustering analysis based on the correlation values at 29 hr. We see that the genes in the selected group show a large fold change and a significant Pvalue in average. For these genes, the enrichment analysis shows that the most significant pathway is influenza A pathway, which shows the involvement of the selected genes in the biological processes of infection.
Figure 4 b presents the dynamical evolution of the gene network respectively at 5, 12, 29, 60 hr for live influenza infection. The selected gene group in Fig. 3 d is placed in the top corner. It can been seen that at 29 hr the structure of the subnetwork of the selected genes changes significantly and thus signals the upcoming deterioration into a disease state which is also in coincidence with the clinic observation.
Therefore, our application results are in coincidence with the experimental observation and successfully detect the earlywarning signal of the impending critical transition.
Discussion and conclusions
Complex diseases significantly damage the health of people all over the world. Detecting the earlywarning signal of the sudden deterioration provides an opportunity to interrupt and prevent the continuing costly cycle of managing these diseases and their complications. Although it is crucial to detect the predisease state so as to prevent the qualitative deterioration by taking appropriate intervention actions, it is a challenging task to reliably identify the predisease state because the state of the system may show neither apparent change nor clear phenomenon before this critical transition during the disease progression. This is also the reason why diagnosis based on traditional biomarkers may fail to indicate a predisease state.
In this work, by detecting the dynamical change of links in a network, we presented a computational method and corresponding algorithm based on HMM to measure the dynamical difference of the system progression, and thus identify the imminent critical transition. It is worth noting that this method aims to detect the earlywarning signal generating from the predisease state (or pretransition state), rather than to find the indication of disease state (or aftertransition state) in which the qualitative deterioration has already taken place.
We applied our method to the identification of the predisease state based on a simulated dataset and two microarray datasets, which demonstrate the sensitivity and effectiveness of our method. For both two diseases, we constructed biomolecular networks (Fig. 4) to gauge the dynamical regulation among genes at different sampling point along a timecourse progression. Both the functional and enrichment analyses validate the computational results. Therefore, the HMMbased method provides a computational possibility of prying into the underlying mechanism of biological processes of the disease progression, and thus may help to achieve the timely intervention. Our dynamic network analysis also suggests, in regard to the diseases, to focus on the specific predisease states to probe the in situ external perturbation (such as environment changes) preceding the development into a badly ill stage. This may lead to not only insights of external environment interactions, but also an effective time window for novel intervention or therapeutic strategies in specific diseases. The main difference between our work and previous ones is that rather than screen out some variables (genes or proteins), the proposed method mainly focuses on the direct identification of critical transition point, by calculating and comparing the consistence probability of each candidate end point of the Markov model in the normal state. Therefore, the accuracy of HMMbased approach is not limited by the selection of features. This is the main value in the potential applications of the HMMbased method from a network point of view.
There are limitations of this work. First, the validity of the identified predisease state and the accurate result needs further supports from animal experiments or clinical studies. Second, the method is insensitive when the correlations are not differentially expressed. Although this work is merely a step towards detecting the earlywarning signals of critical transition during disease progression and the algorithm is expected to be improved in both sensitive and accurate ways, it opens a window of opportunity for the applicable approach to the earlywarning system of critical transition during the biological processes of complex diseases.
Notes
Declarations
Acknowledgements
Publication of this article was funded by National Natural Science Foundation of China (Grant numbers 91530320, 91439103, 11401222 and 61370228); Science and Technology Planning Project of Guangdong Province, China (Grant number 2015B010128008, 2014B090903008, 2015B010109006); Fundamental Research Funds for the Central Universities (Grant number 2014ZZ0064); Pearl River Science and Technology Nova Program of Guangzhou (Grant Number 201610010029).
Declarations
This article has been published as part of BMC Systems Biology Vol 10 Suppl 2 2016: Selected articles from the IEEE International Conference on Bioinformatics and Biomedicine 2015: systems biology. The full contents of the supplement are available online at http://bmcsystbiol.biomedcentral.com/articles/supplements/volume10supplement2.
Author’s contributions
PC and YL conceived the research. PC performed the numerical simulation and real data analysis. All authors wrote the paper and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Venegas JG, Winkler T, Musch G, Melo MF, Layfield D, Tgavalekos N,et al.Selforganized patchiness in asthma as a prelude to catastrophic shifts. Nature. 2005; 434(7034):777–82.View ArticlePubMedGoogle Scholar
 McSharry PE, Smith LA, Tarassenko L.Prediction of epileptic seizures: are nonlinear methods relevantNat Med. 2003; 9(3):241–2.View ArticlePubMedGoogle Scholar
 Pastor BR, Guallar E, Coresh J.Transition models for changepoint estimation in logistic regression. Stat Med. 2003; 22(7):1141–62.View ArticleGoogle Scholar
 Paek SH, Chung HT, Jeong SS, Park C, Kim C, Kim JE, et al.Hearing preservation after gamma knife stereotactic radiosurgery of vestibular schwannoma[J]. Cancer. 2005; 104(3):580–90.View ArticlePubMedGoogle Scholar
 Liu JK, Rovit RL, Couldwell WT. Pituitary apoplexy: Seminars in Neurosurgery. New York: 2001. p. 315–20.Google Scholar
 Tanaka G, Tsumoto K, Tsuji S, Aihara K.Bifurcation analysis on a hybrid systems model of intermittent hormonal therapy for prostate cancer. Physica D: Nonlinear Phenomena. 2008; 237(20):2616–27.View ArticleGoogle Scholar
 Achiron A, Grotto I, Balicer R, Magalashvili D, Feldman A, Gurevich M.Microarray analysis identifies altered regulation of nuclear receptor family members in the predisease state of multiple sclerosis. Neurobiol Dis. 2010; 38(2):201–9.View ArticlePubMedGoogle Scholar
 Chen L, Liu R, Liu ZP, Li M, Aihara K.Detecting earlywarning signals for sudden deterioration of complex diseases by dynamical network biomarkers. Sci Rep. 2012;2. doi:http://dx.doi.org/10.1038/srep00342.
 Liu R, Li M, Liu ZP, Aihara K, Chen L.Identifying critical transitions and their leading biomolecular networks in complex diseases. Sci Rep. 2012;2. doi:http://dx.doi.org/10.1038/srep00813.
 Scheffer M, Bascompte J, Brock WA, Brovkin V, Carpenter SR, Dakos V,et al.Earlywarning signals for critical transitions. Nature. 2009; 461(7260):53–9.View ArticlePubMedGoogle Scholar
 Liu R, Yu X, Liu X, Xu D, Aihara K, Chen L.Identifying critical transitions of complex diseases based on a single sample. Bioinformatics. 2014; 30(11):1579–86.View ArticlePubMedGoogle Scholar
 Litt B, Esteller R, Echauz J, Alessandro MD, Shor R, Henry T, et al.Epileptic seizures may begin hours in advance of clinical onset: a report of five patients. Neuron. 2001; 30(1):51–64.View ArticlePubMedGoogle Scholar
 He D, Liu ZP, Honda M, Kaneko S, Chen L.Coexpression network analysis in chronic hepatitis B and C hepatic lesions reveals distinct patterns of disease progression to hepatocellular carcinoma. J Mol Cell Biol. 2012; 4(3):140–52.View ArticlePubMedGoogle Scholar
 Liu R, Wang X, Aihara K, Chen L.Early diagnosis of complex diseases by molecular biomarkers, network biomarkers, and dynamical network biomarkers. Med Res Rev. 2014; 34(3):455–78.View ArticlePubMedGoogle Scholar
 Liu R, Aihara K, Chen L.Dynamical network biomarkers for identifying critical transitions and their driving networks of biologic processes. Quant Biol. 2013; 1(2):105–14.View ArticleGoogle Scholar
 Liu X, Liu R, Zhao XM, Chen L.Detecting earlywarning signals of type 1 diabetes and its leading biomolecular networks by dynamical network biomarkers. BMC Med Genom. 2013; 6(Suppl 2):S8.Google Scholar
 Tan Z, Liu R, Zheng L, Hao S, Fu C, Li Z, et al.Cerebrospinal fluid protein dynamic driver network: At the crossroads of brain tumorigenesis. Methods. 2015; 83:36–43.View ArticlePubMedGoogle Scholar
 Zeng T, Zhang C, Zhang W, Liu R, Liu J, Chen L.Deciphering early development of complex diseases by progressive module network. Methods. 2014; 67(3):334–43.View ArticlePubMedGoogle Scholar
 Liu R, Chen P, Aihara K, Chen L.Identifying earlywarning signals of critical transitions with strong noise by dynamical network markers. Sci Rep. 2015;5. doi:http://dx.doi.org/10.1038/srep17501.
 Li M, Zeng T, Liu R, Chen L.Detecting tissuespecific early warning signals for complex diseases based on dynamical network biomarkers: study of type 2 diabetes by crosstissue analysis. Brief Bioinformatics. 2014; 15(2):229–43.View ArticlePubMedGoogle Scholar
 Chen P, Liu R, Chen L, Aihara K.Identifying critical differentiation state of MCF7 cells for breast cancer by dynamical network biomarkers. Front Genet. 2015;6. doi:http://dx.doi.org/10.3389/fgene.2015.00252.
 Chen L, Wang RS, Zhang XS. Biomolecular networks: methods and applications in systems biology. New Jersey: John Wiley & Sons; 2009.View ArticleGoogle Scholar
 Wurmbach E, Chen Y, Khitrov G, Zhang W, Roayaie S, Schwartz M, et al.Genomewide molecular profiles of HCVinduced dysplasia and hepatocellular carcinoma. Hepatology. 2007; 45(4):938–47.View ArticlePubMedGoogle Scholar
 Bruix J, Boix L, Sala M, Llovet JM. Focus on hepatocellular carcinoma. Cancer Cell. 2004; 5(3):215–9.View ArticlePubMedGoogle Scholar
 Farazi PA, DePinho RA. Hepatocellular carcinoma pathogenesis: from genes to environment. Nat Rev Cancer. 2006; 6(9):674–87.View ArticlePubMedGoogle Scholar
 Huang Y, Zaas AK, Rao A, Dobigeon N, Woolf PJ, Veldman T, et al.Temporal dynamics of host molecular responses differentiate symptomatic and asymptomatic influenza a infection. PLoS Genet. 2011; 7(8):e1002234.View ArticlePubMedPubMed CentralGoogle Scholar