Identifying the critical state of cancers by single-sample Markov flow entropy

Theoretical background

The DNB theory (Chen et al., 2012; Liu, Aihara & Chen, 2013) is proposed to visualize the evolution of complex diseases into three periods: a normal (before-transition) state, a pre-disease (critical) state, and a disease (after-transition) state. The normal state is stable and resilient to perturbations, while the pre-disease state is unstable and reversible with appropriate intervention. The disease state is also stable and resilient, but is almost irreversible (Litt et al., 2001; Liu, Chen & Chen, 2020; McSharry, Smith & Tarassenko, 2003; Venegas et al., 2005). Moreover, when a complex system is near the critical point, a dominant group of DNB biomolecules can be identified, based on observed data, that satisfy the following three conditions (Chen et al., 2012a):

• The Pearson correlation coefficient (PCC_in) between any pair of members in the DNB group rapidly increases;

• The Pearson correlation coefficient (PCC_out) between one member of the DNB group and any other non-DNB member rapidly decreases;

• The standard deviation (SD_in) or coefficient of variation for any member in the DNB group drastically increases.

Our DNB theory suggests that phase transitions are characterized by a set of highly fluctuating and highly correlated features/variables, which implies an imminent transition to the disease state. Therefore, we can quantify the tipping points based on these three conditions, provide early-warning signals of disease deterioration, and then determine the biomolecular dominant group to make up DNB members. These three conditions are the basis of DNB theory and have many applications in disease progression and biological processes to predict its critical state (Liu, Aihara & Chen, 2013; Liu, Chen & Chen, 2020).

Markov chain entropy (MCE) (Shi et al., 2018) was constructed to represent the potential of samples containing scRNA-Seq or bulk RNA-Seq information. This method requires the normalized RNA Seq profile and a defined interaction network between the genes.

Specifically, given an interaction network G with $n_v$ nodes and $n_E$ edges, and the transition probability matrix $P=\left\{p_{ij}\right\}$, where $i,j=1,\dots ,n_E$. Where $p_{ij}=0$ if there is no edge between two different nodes i and j. Then the Markov chain entropy based on network G is defined as

(1) $$\mathrm{M}\mathrm{C}\mathrm{E}=-\sum _{(\mathrm{i},j)\in \overline{E}}{\pi }_i{p}_{ij}\log \left({\pi }_i{p}_{ij}\right)$$where ${\pi }_i$ was node i’s expression normalized to [0,1]. In this article, MCE was split into each node to form Markov flow entropy (MFE) (Fig. 1B). For centered gene $g^k$, the MFE is defined as

(2) $$\mathrm{M}\mathrm{F}\mathrm{E}(g^k)=-\sum _{j=1}^M{\pi }_i{p}_{ij}\log \left({\pi }_i{p}_{ij}\right)$$where M denotes the number of the 1st-order neighbors ${(\mathrm{g}}_1^k,\dots ,{\mathrm{g}}_M^k)$ of $g^k$, $p_{ij}$ denotes the transition probability from $g^k$ to ${\mathrm{g}}_j^k$, and ${\pi }_i{p}_{ij}$ means the probability transmitted from $g^k$ to ${\mathrm{g}}_j^k$ along the directed edge.

Algorithm to identify the tipping point based on sMFE

Reference samples (healthy or relatively healthy cells; tumor-adjacent samples here) and only one disease sample were required in the sMFE method to identify the tipping point with the following algorithm (Fig. 1B).

(Step 1) Constructing a global directed network ${\mathrm{N}}^G$ by mapping the genes to a protein–protein interaction (PPI) network from STRING (http://string-db.org) (Szklarczyk et al., 2015), with a confidence level of 0.800. Discard isolated nodes that do not have any links with other nodes. The direct interaction network was constructed based on rewiring the protein–protein interaction (PPI) network by GNIPLR (gene networks inference using projection and lagged regression) (Zhang, Chang & Liu, 2021) method. Tumor-adjacent samples were used to infer the gene regulatory network.

(Step 2) Mapping data into the global network ${\mathrm{N}}^G$ constructed in the previous step to form expression matrix. Mapping datasets of four cancers downloaded from the TCGA database to the global network ${\mathrm{N}}^G$, then the gene expression matrix formed for further analysis.

(Step 3) Calculating the local Markov flow entropy (MFE) for each gene. Specifically, we extract each gene $g^k$’s local network ${\mathrm{N}}^k$ (k = 1, 2, …, L) from the global network ${\mathrm{N}}^G$ (Fig. 1B), the number of the local network is L, and $\left\{{\mathrm{g}}_1^k,\dots ,{\mathrm{g}}_M^k\right\}$ are the 1st-order neighbors of $g^k$. Then, based on n reference samples, the formula for calculating local Markov flow entropy $MF{E}_k^n(t)$ is

(3) $$MF{E}_k^n(t)=-\sum _{i=1}^M\overline{x_k^n}{p}_{ki}^n\log \left(\overline{x_k^n}{p}_{ki}^n\right),$$with

(4) $$p_{ki}^n(t)={\displaystyle \frac{\left|\mathrm{P}\mathrm{C}{\mathrm{C}}^n\left(g_i^k(t),g^k(t)\right)\right|}{\sum _{j=0}^M\left|\mathrm{P}\mathrm{C}{\mathrm{C}}^n\left(g_j^k(t),g^k(t)\right)\right|},}$$where the constant M represents the count of neighbors in the local network ${\mathrm{N}}^k$ and $PC{C}^n\left(g_i^k(t),g^k(t)\right)$ represents the Pearson correlation coefficient (PCC) between the center gene ${\mathrm{g}}^k$ and its 1-st neighbor ${\mathrm{g}}_i^k$ based on n reference samples at time point t. Specifically, ${\mathrm{g}}_i^k{=\mathrm{g}}_{}^k$ if $j=0$. $\overline{x_k^n}$ is the mean of the canter gene $g_{}^k$’s normalized expression based on n reference samples at time point t.

(Step 4) Calculating the differential MFE called $sMF{E}_k^{}(t)$. Under n + 1 mixed samples, where a single sample from an individual was added to n reference samples at time point t, the same with the prior step, $MF{E}_k^{n+1}(t)$ is calculated, i.e.,

(5) $$MF{E}_k^{n+1}(t)=-\sum _{i=1}^M\overline{x_k^{n+1}}{p}_{ki}^{n+1}\log \left(\overline{x_k^{n+1}}{p}_{ki}^{n+1}\right)$$

Then, the $sMF{E}_k^{}(t)$ is calculated as

(6) $$sMF{E}_k(t)={\displaystyle \frac{1}{\log M}\left|S{D}^n\left(g^k(t)\right)\ast MF{E}_k^{nt}(t)-S{D}^{n+1}\left(g^k(t)\right)+MF{E}_k^{n+1}(t)\right|,}$$where $\mathrm{S}{\mathrm{D}}^n(g^k(t))$ and $\mathrm{S}{\mathrm{D}}^{n+1}(g^k(t))$ represent the standard deviation of the gene expression for the center gene ${\mathrm{g}}^k$ based on n reference samples and n + 1 mixed samples at time point t, respectively.

(Step 5) The fifth step is to calculate the global single-sample Markov flow entropy $sMFE(t)$, i.e.,

(7) $$sMFE(t)={\displaystyle \frac{1}{Q}\sum _{k=1}^{Q}sMF{E}_k(t),}$$where constant $Q$ represents the count of top 5% genes with highest local sMFE. Furthermore, t-test was applied to test if the identified critical state significantly different with the other states. The alternative state was viewed as critical state if P < 0.05 (see Supplemental Information). From the formula of sMFE, we can see that when the interaction between the center node and the 1-st neighbor node converges (that is, the transition probability is the same), the sMFE value reaches the maximum value. Conversely, when the interaction between the center node and one 1-st neighbor node is extremely strong (the transition probability equals 1), the sMFE value reaches the minimum. From DNB theory, it can know that in the sample at the critical point period, the interaction between DNB molecules and their neighboring nodes will increase extremely so that it can form a significant difference from the reference sample, and sMFE can measure this difference, and when the sample reaches the critical point, the significance increases so that it can warn the disease evolution process. According to DNB theory (Chen et al., 2012; Liu, Aihara & Chen, 2013), when the system approaches the critical point, the sMFE score can effectively characterize the fluctuation of the network, and thus can be regarded as an early warning signal prior to critical deterioration.

Data processing and functional analysis

The Cancer Genome Atlas (TCGA) database (GDC (cancer.gov)) was used to download four unrelated clinical tumor datasets: esophageal carcinoma (ESCA), colon adenocarcinoma (COAD), (kidney renal clear cell carcinoma (KIRC)), and lung adenocarcinoma (LUAD). RNA-seq data from tumor and tumor-adjacent samples, as well as clinical information, were included in these datasets. Then we classified the tumor samples into different stages based on clinical (stage) information downloaded from TCGA, samples lacking stage information were discarded, and Table S1 provided the full clinical staging information.

The molecular global template network was constructed using the steps shown in Fig. 1A.

The first step in constructing the molecular global template network was to download the protein–protein interaction networks for Homo sapiens from STRING (http://string-db.org). This data was then used to generate a global template network, which was composed of nodes representing proteins and edges representing interactions between proteins. All interaction used for discussion was picked with high confidence (higher than 0.8). Nodes that are not connected to any other nodes are removed from the PPI network.

Second, the direct interaction network was constructed based on rewiring the protein–protein interaction (PPI) network by GNIPLR (gene networks inference using projection and lagged regression) (Zhang, Chang & Liu, 2021) method. Tumor-adjacent samples were used to infer the gene regulatory network (GRN; directed network) for subsequent analysis (detailed information for how to choose adjustable parameters in GNIPLR were added in Supplemental Information).

In the end, the directed network, created in prior step, was adapted to map the gene expression of each RNA-seq dataset to form a molecular network for further analysis.

The Kyoto Encyclopedia of Genes and Genomes (KEGG) database Mapper tool (KEGG Mapper Color) and the Database for Annotation, Visualization and Integrated Discovery (DAVID) Functional Annotation Tool (DAVID: Functional Annotation Tools (ncifcrf.gov)) were used to conduct the enrichment analysis. These statistical analysis carried in the article was done using R Software v 4.0.3 (R: The R Project for Statistical Computing (r-project.org)), and the R Survival package was used to make Cox survival analysis. The visualization results and network analysis are implemented using Cytoscape software (www.cytoscape.org).

Identifying the critical transition points of four cancers with sMFE

The sMFE algorithm was used to identify the critical transition points of four cancers (ESCA, COAD, KIRC, and LUAD) from TCGA datasets. These datasets contain much clinical information, in which doctors divide cancer patients into different stages according to the size of the tumor diameter, such as KIRC samples were divided into four stages (I, II, III, and IV), while LUAD samples were divided into seven stages (IA, IB, IIA, IIB, IIIA, IIIB, and IV). The tumor-adjacent (TA) samples were deemed relatively normal and used as a reference for each cancer when calculating sMFE according to the proposed algorithm. The average sMFE score curve for each cancer stage was calculated to reflect the dynamic change and was shown in Figs. 2A–2D. Based on DNB theory, the critical state during disease progression can be reflected by the tipping point of the sMFE score curve (Fig. 1A). Kaplan–Meier (log-rank) survival analysis was performed on the prognosis of before-transition and after-transition samples (Figs. 2E–2H) to further validate the identified critical state.

As seen in Fig. 2A (yellow curve), the average sMFE score of ESCA increased significantly (P = 0.0077) at stage IIIA, suggesting an upcoming critical transition after stage IIIA.; that is, the tumor migrates and invades to form distant metastasis at stage IIIB and ultimately cause distant metastasis (stage IV) (Dong et al., 2022). While at the six-time points, there was little significant difference among the average gene expression of differentially expressed genes (DEGs) (the gray curve in Fig. 2A). Of note, samples taken before and after the transition point (Stage IIIA) had significantly different survival curves (P < 0.0001; Fig. 2E), and patients taken after the transition point (Stage IIIA) had shorter survival times than those before the transition.

When applied to COAD, a significant change (P = 0.0045; Fig. 2B) in sMFE was detected around stage IIB, suggesting lymph node metastasis and tumor invasion of other adjacent organs in stage III (Hari et al., 2013). As presented in gray in Fig. 2B, dynamic changes (P = 0.16) in the mean gene expression of DEGs were detected in stage IIIB and therefore failed to provide timely early-warning signals of the critical transition. Moreover, it is worth noting that survival analysis showed that the survival time of samples taken before the transition point (Stage IIB) was significantly different than that taken after the transition point (Stage IIB) (P < 0.0001; Fig. 2F). Furthermore, patients from after the identified critical state (Stage IIB) had shorter survival times than those from before the transition.

For KIRC dataset, a drastic increase (P = 0.014) in the sMFE score was observed in stage II (transition point; Fig. 2C), indicating an upcoming critical transition into the disease state (stage III/IV), where the tumor migrates to form distant metastasis at stage III (Gu & Zhao, 2019). In terms of mean gene expression, dynamic changes (P = 0.25; the gray curve in Fig. 2C) were detected in stage III and therefore failed to provide timely early-warning signals of the critical transition. Furthermore, survival analysis demonstrated that survival curves had a significant difference (P < 0.0001) between stage I~II and stage III~IV KIRC samples according to clinical information (Fig. 2G). In addition, the survival time of patients after stage II was significantly shorter than before stage II.

The red curve presented the LUAD data in Fig. 2D, the drastic increase of sMFE occurred at stage IIIB (P = 0.0027), indicating an upcoming critical transition after stage IIIB; that is, stage IV was characterized by a distant metastasis process, in which the tumor cells invaded distant tissues or organs (Chiang & Massagué, 2008). However, the gray curve in Fig. 2D shows little significant increase in the mean expression of DEGs occurring among the six-time points. Of note, the survival curves of samples taken before and after the transition point (stage IIIB) tend to show a significantly different (P = 0.0024; Fig. 2H). In addition, patients with stage after stage IIIB had shorter survival times than those before. As shown, sMFE is an effective approach to identify the critical state of cancers.

The dynamic evolution of gene regulatory networks

To validate the critical states identified by sMFE, we further examined them from the network perspective. First, for each sample, the top 5% of genes with the highest sMFE scores were considered the signaling genes. Dynamic network biomarkers (DNBs) consist of common signaling genes in samples in the identified critical state. In DNB theory, these DNB molecules exhibit statistical properties when the system reaches a critical point. Biologically, they may be involved in important biological reaction processes during disease progression. These DNB molecules were mapped to the STRING database to form a connected PPI network, and their dynamic changes in the process of disease evolution were studied from the network level.

The landscape of the local sMFE score for ESCA data was depicted in Fig. 3A. It was worth noting that the local sMFE scores of the signaling genes showed an abruptly increased collaborative manner around stage IIIA. Furthermore, in the network structure, a noticeable change occurred at stage IIIA (Fig. 3E), indicating an impending critical transition (Dong et al., 2022), which spontaneously come together with the experimental results. When applied to the COAD dataset, the landscape of the local sMFE score was shown in Fig. 3B. It was discovered that the peak of local sMFE scores for signaling genes appeared at stage IIB. Furthermore, there was a significant change in the network structure at stage IIB (Fig. 3F), indicating lymph node metastasis and tumor invasion of other adjacent organs in stage III (Hari et al., 2013). As shown in Fig. 3C, for KIRC, showing that the peak local sMFE of signaling genes appeared at stage II. Moreover, in the PPI network, a drastic change occurred in stage II (Fig. 3G), where the tumor gradually migrates to form distant metastasis at stage III (Gu & Zhao, 2019). Results for LUAD in Fig. 3D, the local sMFE score landscape peaked in stage IIIB and there was a significant change in the network structure of stage IIIB (Fig. 3H), i.e., stage IV was characterized by a distant metastasis process, in which the tumor cells invaded distant tissues or organs (Chiang & Massagué, 2008). In addition, Fig. S2 provided results of the dynamic network changes across all stages.

Prognostic biomarkers of tumors in the sMFE method

In order to identify biomarkers that were effective in predicting prognoses. Indeed, signaling genes could be categorized into two types of molecules for prognostic prediction as common biomarkers for all samples: those for samples with poor prognosis, termed pessimistic sMFE (P-sMFE) biomarkers; and those with good prognosis, termed optimistic sMFE (O-sMFE) biomarkers. Furthermore, those two types of biomarkers are not only the signaling genes, pessimistic/optimistic biomarkers, but also the non-differentially expressed genes usually screened out by most researchers in the first step. More details about the identification method are provided in Fig. S1.

The meaning of these two biomarkers can also be understood in this way, samples of optimistic sMFE biomarkers that appeared in their signaling gene had better-predicted prognoses than those without; in other words, these samples were expected to survive longer than others. Conversely, samples with pessimistic sMFE biomarkers in their signaling genes had lower predicted prognoses than those without; that is, these samples were expected to survive for a shorter time than other samples. A total of four genes were identified as optimistic sMFE biomarkers for ESCA (Fig. 4A and Table S2), while six genes were identified as pessimistic sMFE biomarkers (Fig. 4A and Table S2). Specifically, survival times of samples with O-sMFE biomarkers SCYL3 and TRPS1 (P = 0.016 and 0.036 respectively) were longer than those without, as shown in Fig. 4A. Conversely, samples with P-sMFE biomarkers EPCAM and POLK (P = 0.0057 and 0.0064, respectively) had shorter survival times than those without.

In the analysis for COAD, AGR2 (P = 0.044; Fig. 4B) and PCM1 (P = 0.046; Fig. 4B) were identified as O-SMFE biomarkers, and it’s seen the survival time of these identified samples tend to be longer than that of other samples according to survival curve analysis (Fig. 4B). Additionally, ATP1A3 (P < 0.0001; Fig. 4B) and PARL (P = 0.00093; Fig. 4B) were identified as P-SMFE biomarkers, survival analysis in Fig. 4B showed the survival time of these identified samples tend to be shorter than that of other samples. When applied to KIRC dataset, samples with O-sMFE biomarkers MAML2 and RNF5 (P = 0.0058 and 0.0034, respectively) had better prognoses, with longer survival times, compared to other samples. Conversely, samples with P-sMFE biomarkers BST1 and VPS39 (P = 0.0085 and 0.00011, respectively) had worse prognoses, with shorter median survival times of almost 5 years. As shown in Fig. 4D, for LUAD, ACOX2 (P = 0.021; Fig. 4D) and KHDRBS3 (P = 0.027; Fig. 4D) were identified as O-SMFE biomarkers and FLNC (P = 0.0021; Fig. 4D) and FZD1 (P = 0.0091; Fig. 4D) were identified as P-SMFE biomarkers.

We propose a method to uncover novel biomarkers, drug targets and key regulators, which are usually overlooked by traditional studies due to their non-differential expression. In addition, all identified optimistic/pessimistic sMFE biomarkers for the four cancers were supplemented in Table S2.

Functional analysis of common signaling genes

Common signaling genes for the four cancers were extracted for functional analysis (Fig. 5A). As shown in Fig. 5B, not only were there occurred rich intersections among signaling genes from different cancers, but their biological functions were also closely related. In addition, we performed the KEGG pathway enrichment analysis of common signaling genes, where their samples located in tipping points, to explore the biological processes in which they participated. For these four tumor datasets, as shown in Figs. 5C–5F, these signaling genes were significantly associated with cancer-associated pathways. For example the cell cycle (Hartwell & Kastan, 1994), PI3K-Akt signaling pathway (Martini et al., 2014), spliceosome (Yang, Beutler & Zhang, 2022), Wnt signaling pathway (Zhou et al., 2022). In addition, Figs. 5G and 5H illustrated that common signaling genes were mainly enriched in the above-mentioned pathway, which suggests that these common signaling genes may serve an important contributor to cancer development and progression. For instance, Ras-related C3 botulinum toxin substrate 2 (RAC2) promotes abnormal proliferation of quiescent cells by enhanced JUNB expression via the MAL-SRF pathway (Pei et al., 2018), PIK3CA-mutation was reported related to metastatic breast cancer (Higgins et al., 2012; Mosele et al., 2020), Protein tyrosine kinase 2 (PTK2) was a novel therapeutic target to overcome acquired EGFR-TKI resistance in non-small cell lung cancer and as a driver of radioresistance in HPV-negative head and neck cancerPTK2/FAK and radioresistance (Skinner et al., 2016; Tong et al., 2019), RhoA is associated with invasion and poor prognosis in colorectal cancer (Jeong et al., 2016).