CDNFE suggests FNDC3B and NECTIN4 as drivers of precancer progression via PI3K/AKT EMT.

The CDNFE method was applied to two real datasets. The first dataset from Luohe Central Hospital contains a total of 9 samples with uterine prolapse or endometriosis hysterectomy and immunosequencing data. To construct a high-quality single-cell GEM for subsequent tipping point detection, we performed systematic filtering of the raw sequencing data following widely adopted quality control practices in the field. Stringent threshold-based filtering was applied to remove low-quality cells and potential doublets: cells with either too few (5000) detected genes were excluded, and those with mitochondrial gene content exceeding 20% were identified as compromised or dying cells and removed. Furthermore, genes detected in extremely few cells were filtered out, retaining only those expressed in at least three cells to ensure analytical robustness. All thresholds were tailored according to the specific experimental system and data distribution to preserve biologically relevant cell populations while minimizing technical noise. Among the nine samples, three were identified as HPV infections, three as precancerous lesions, and three as cervical cancer. In our study, HPV infections were classified as Normal (3 samples), high-grade squamous intraepithelial lesions as HSIL (3 samples), and cervical cancer was further divided into early CESC (2 samples) and late CESC (1 sample). For single-cell sequencing data, these samples were subjected to 10X single-cell RNA-seq sequencing, data quality control, cellranger quantitative analysis. Firstly, the samples were subjected to cell quality control, and the cell concentration was adjusted to the ideal concentration for subsequent sequencing. The second step was the preparation of cDNA sequencing libraries and the construction of cDNA libraries. The third step is high-throughput sequencing on the machine. Finally, the raw data obtained will be split according to the unique library structure of 10xscRNA-Seq, and the Barcode, UMI, and insertion fragment parts of the reads will be split, after which the insertion fragment parts will be compared to the reference genome, and then the ratio of the comparison to each region will be counted and the expression amount will be counted. After filtering, 75758 cells (8077/7775/5366 for normal samples; 13256/9466/9583 for HSIL samples; 7018/7211/8006 for CESC samples) were retained for the following analysis. The human data used in this study were obtained from the Luohe Central Hospital and were approved by the Ethics Committee of Luohe Central Hospital (Approval No. Z20221343023). The study was conducted in accordance with the Declaration of Helsinki and relevant national guidelines for biomedical research involving human participants. Written informed consent was obtained from all participants for the use of their samples and associated data for research purposes. We analyzed the preprocessed gene expression datasets measured in CESC. The second dataset was downloaded from the NCBI GEO database (access ID: GSE63514 ). The dataset comprised 54675 probes. Then, we transformed the downloaded matrixes into required GEMs through ID conversion and the deletion of duplicate genes and null values. In the case of multiple probes corresponding to the same gene, the values were individually averaged for each to obtain GEMs. The dataset was composed of 128 samples, including 24 normal specimens, 14 CIN1, 22 CIN2, 40 CIN3 and 28 cancer specimens, where CIN1, 2, and 3 corresponded to low-grade lesions, high-grade lesions, and carcinoma in situ, respectively. We classified 24 normal samples as Normal, 14 CIN1 samples as low-grade squamous intraepithelial lesion (LSIL), and grouped 22 CIN2 and 40 CIN3 samples as HSIL according to cervical cancer standards. Additionally, 28 cancer samples were categorized as CESC. For pseudotime analysis, trajectory inference was performed using the Monocle package in R. A CellDataSet object was constructed from the normalized expression data, specifying negbinomial.size() as the expression family. Ordering genes was identified through differential gene expression testing using the model formula ~ celltype. Dimensionality reduction was carried out using the DDRTree method with max_components = 2 and num_dim = 10. Pseudotime was subsequently calculated using the orderCells function without manual specification of the root cell; instead, the cell state with the lowest pseudotime value was automatically selected as the root by the algorithm. For ROC analysis, binary labels were assigned to samples based on whether they belonged to the critical state (HSIL). ROC curves were generated using the pROC R package, and AUC values were calculated via resampling. The enrichment analysis of DNBs is based on the Gene Ontology Consortium ( http://geneontol-ogy.org ), DAVID Bioinformatics Resources ( https://david.ncifcrf.gov/ ) and Circos ( http://www.circos.ca/ ). The gene function annotation of each disease was obtained through GeneCards ( http://www.genecards.org/ ). Protein–Protein Interaction (PPI) networks were drawn using STRING ( https://string-db.org/ ) and the client software Cytoscape ( https://cytoscape.org/ ). The public dataset 52 comprises 38 samples, including 10 normal (HPV-negative), 7 HSIL, and 21 CESC. Our dataset consists of 9 samples, including 3 normal (HPV-positive), 3 HSIL, and 3 CESC. To integrate the two datasets, we first extracted the overlapping gene set and then applied the ComBat algorithm to correct for batch effects. After merging, the combined cohort comprised four groups: 10 normal (HPV-negative), 3 normal (HPV-positive, LSIL), 10 HSIL, and 24 CESC. We then applied the CDNFE model to the expanded cohort to identify the critical transition point. The results confirmed that the critical point remains localized to the HSIL stage, consistent with findings from the original cohort alone. The CDNFE method is fundamentally designed to detect critical or pre-disease states in the progression of cervical cancer, with a particular focus on identifying DNBs. The schematic representation of the CDNFE algorithm is depicted in Fig. 1 , providing a visual overview of the process. This approach will be explored in this section, detailing each step to provide a comprehensive understanding of the CDNFE approach. [Step 1] Construct the directed network at the specified time point t . Utilizing the WGCNA co-expression network along with gene expression data, a directed network was meticulously constructed based on the directional index. This index was delineated as follows: 1 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\omega }_{i,j}=\mathop{\sum }\limits_{\hat{x}\in \vec{\hat{X}}}\mathop{\sum }\limits_{y\in \vec{Y}}p\left(\hat{x},y\right)log\frac{p\left(\hat{x},y\right)}{p\left(\hat{x}\right)p\left(y\right)}-\mathop{\sum }\limits_{\hat{x}\in \vec{X}}\mathop{\sum }\limits_{y\in \vec{Y}}p\left(x,y\right)log\frac{p\left(x,y\right)}{p\left(x\right)p\left(y\right)}$$\end{document} ω i , j = ∑ x ˆ ∈ X ˆ ⃗ ∑ y ∈ Y ⃗ p x ˆ , y l o g p x ˆ , y p x ˆ p y − ∑ x ˆ ∈ X ⃗ ∑ y ∈ Y ⃗ p x , y l o g p x , y p x p y Vectors \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vec{U}$$\end{document} U ⃗ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vec{V}$$\end{document} V ⃗ represent the expression profiles of genes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}_{i}$$\end{document} g i and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}_{j}$$\end{document} g j in the sample/cell, respectively, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vec{\hat{X}}$$\end{document} X ˆ ⃗ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vec{X}$$\end{document} X ⃗ are defined as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vec{\hat{X}}=\left(\vec{U}+\vec{V}\right)/\sqrt{2}$$\end{document} X ˆ ⃗ = U ⃗ + V ⃗ / 2 , which ensures that information from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,\vec{U}$$\end{document} U ⃗ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vec{V}$$\end{document} V ⃗ has equal weight and prevents the exponential growth of the feature values during the crossover process and improves the numerical stability of the algorithm, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vec{X}=\vec{U}$$\end{document} X ⃗ = U ⃗ , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{Y}$$\end{document} Y ˆ is the phenotype representing the binary vector for each sample (0-1). \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(\hat{x},y)$$\end{document} p ( x ˆ , y ) represents the joint probability density function (pdf) of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vec{\hat{X}}$$\end{document} X ˆ ⃗ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vec{Y}$$\end{document} Y ⃗ , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(x,y)$$\end{document} p ( x , y ) represents the pdf of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vec{X}$$\end{document} X ⃗ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vec{Y}$$\end{document} Y ⃗ , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p\left(\hat{x}\right)$$\end{document} p x ˆ , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p\left(x\right)$$\end{document} p x , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(y)$$\end{document} p ( y ) represent the edge pdf of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vec{\hat{X}}$$\end{document} X ˆ ⃗ , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vec{X}$$\end{document} X ⃗ , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vec{Y}$$\end{document} Y ⃗ , respectively. The positive determination value indicates that the integration of the gene \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}_{j}$$\end{document} g j is an improvement of the gene \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}_{i}$$\end{document} g i mutual information (MI), that is, in the directional network, there is a directed edge \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({g}_{i},{g}_{j})$$\end{document} ( g i , g j ) from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}_{i}$$\end{document} g i to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}_{j}$$\end{document} g j . When the direction index is greater than 0, we can determine the existence of an edge from one gene to another. The theoretical underpinnings of network construction are detailed in the Files S2 and S3 . [Step 2] Extract the local network of each gene. In our approach, each local network considers exclusively the first-order and second-order neighbors. Local network \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${M}^{k}$$\end{document} M k centers on gene \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}^{k}$$\end{document} g k , which has a first-order in-degree neighbors \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{{g}_{i,1}^{k},{g}_{i,\,2}^{k},\cdots ,{g}_{i,b}^{k}\}$$\end{document} { g i , 1 k , g i , 2 k , ⋯ , g i , b k } and b first-order out-degree neighbors \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{{g}_{o,1}^{k},{g}_{o,2}^{k},\cdots ,{g}_{o,b}^{k}\}$$\end{document} { g o , 1 k , g o , 2 k , ⋯ , g o , b k } . [Step 3] Filtering the local directed network using cross-over strength. Subsequently, we filter our extracted local directed networks based on the cross-over strength, and the formula for calculating the crossover strengths for the first-order neighborhood gene \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{g}_{i}}^{k}$$\end{document} g i k is defined as follows: 2 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,{D}_{i}=\lambda {D}_{i}^{{in}}+\left(1-\lambda \right){D}_{i}^{{out}}$$\end{document} D i = λ D i in + 1 − λ D i out where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document} λ is a constant, 0 <  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document} λ  < 1, and the value can be ascertained on an individual, case-specific basis. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${D}_{i}^{{in}}={\sum }_{j=1}^{{a}_{i}}{\omega }_{{ij}}$$\end{document} D i in = ∑ j = 1 a i ω ij is the in-strength of the first-order neighborhood genes, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${a}_{i}$$\end{document} a i represents the number of second-order in-degree neighbors, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${D}_{i}^{{out}}={\sum }_{j=1}^{{b}_{i}}{\omega }_{{ij}}$$\end{document} D i out = ∑ j = 1 b i ω ij is the out-strength of the first-order neighborhood genes, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${b}_{i}$$\end{document} b i represents the number of second-order out-degree neighbors. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\omega }_{{ij}}$$\end{document} ω ij represents the weight of the directed edge \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({g}_{i},{g}_{j})$$\end{document} ( g i , g j ) , which is defined as 3 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\omega }_{{ij}}=\frac{{W}_{{ij}}}{{\sum }_{j=1}^{{a}_{i}/{b}_{j}}{W}_{{ij}}},{W}_{{ij}}=\frac{\left(n-1\right)\varDelta {pcc}\left({g}_{i}^{k},{g}_{j}^{k}\right)}{1-{\left({{pcc}}^{n}\left({g}_{i}^{k},{g}_{j}^{k}\right)\right)}^{2}},$$\end{document} ω ij = W ij ∑ j = 1 a i / b j W ij , W ij = n − 1 Δ pcc g i k , g j k 1 − pcc n g i k , g j k 2 , based on our calculated crossover strengths, the first-order neighborhood genes are ranked, and the first-order neighborhood genes ranked in the top N are finally selected, and the local network containing the first-order neighborhood genes is saved. [Step 4] Determine the local causality directed network. Using the local directed network obtained in step 3, we construct the causal directed network of the samples based on the causal strength, if the causal strength is greater than 0, then there are causal directed edges; otherwise, there are none. The formula for the causal intensity metric is as follows: 4 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${w}_{t}\left({g}^{k},{g}_{j}^{k}\right)={\mathrm{ln}}\left(\frac{\hat{\beta }}{\beta }\right),$$\end{document} w t g k , g j k = ln β ˆ β , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\beta }$$\end{document} β ˆ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} β denote the test errors obtained from Eqs. 5 and 6 , respectively, when applied to the case sample. 5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${x}^{k}=\hat{f}\left({x}_{1}^{k},{x}_{2}^{k},\cdots ,{x}_{j-1}^{k},{x}_{j+1}^{k},\cdots ,{x}_{N}^{k}\right)+\hat{\beta },$$\end{document} x k = f ˆ x 1 k , x 2 k , ⋯ , x j − 1 k , x j + 1 k , ⋯ , x N k + β ˆ , 6 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${x}^{k}=f\left({x}_{1}^{k},{x}_{2}^{k},\cdots ,{x}_{j-1}^{k},{x}_{j}^{k},{x}_{j+1}^{k},\cdots ,{x}_{N}^{k}\right)+\beta$$\end{document} x k = f x 1 k , x 2 k , ⋯ , x j − 1 k , x j k , x j + 1 k , ⋯ , x N k + β where symbol \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${x}_{j}^{k}$$\end{document} x j k represents the expression of the gene \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}_{j}^{k}$$\end{document} g j k of the local directed network \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${M}^{k}$$\end{document} M k centered with the gene \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}^{k}$$\end{document} g k . Specifically, for a local network \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${M}^{k}$$\end{document} M k centered with a gene \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}^{k}$$\end{document} g k in the network, whose first-order neighbors are genes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\{{g}_{1}^{k},{g}_{2}^{k},\cdots ,{g}_{N}^{k}\right\}$$\end{document} g 1 k , g 2 k , ⋯ , g N k , we assume that the change of expressions of any first-order neighbor may affect that of the center node \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}^{k}$$\end{document} g k . A group of relative healthy reference samples/cells serves as the training samples/cell to determine regression models \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{f}$$\end{document} f ˆ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f$$\end{document} f , while a single case sample/cell at each time point t is designated as the test sample/cell. By inputting the test sample, into \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{f}$$\end{document} f ˆ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${f}$$\end{document} f . [Step 5] Calculation of a local CDNFE score for each local causality directed network. We can obtain the local causality directed network \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${M}^{k}\left(k=\mathrm{1,2},\cdots ,m\right)$$\end{document} M k k = 1,2 , ⋯ , m , and each causality local directed network is centered at gene \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}^{k}$$\end{document} g k with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L$$\end{document} L first-order out-degree neighbors \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{{g}_{o,1}^{k},{g}_{o,2}^{k},\cdots ,{g}_{o,L}^{k}\}$$\end{document} { g o , 1 k , g o , 2 k , ⋯ , g o , L k } and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W$$\end{document} W first-order in-degree neighbors \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{{g}_{i,1}^{k},{g}_{i,2}^{k},\cdots ,{g}_{i,W}^{k}\}$$\end{document} { g i , 1 k , g i , 2 k , ⋯ , g i , W k } , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${W}+L=N$$\end{document} W + L = N , its corresponding local CDNFE score at t is defined below. 7 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,{\,{CDNFE}}_{t}^{K}=\frac{W}{N}{{CDNFE}}_{i}^{k}+\frac{L}{N}{{CDNFE}}_{o}^{k},$$\end{document} CDNFE t K = W N CDNFE i k + L N CDNFE o k , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{CDNFE}}_{{in}}^{k}$$\end{document} CDNFE in k and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{CDNFE}}_{{out}}^{k}$$\end{document} CDNFE out k are defined as follows, 8 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,{{CDNFE}}_{i}^{k}=-\frac{1}{n}\frac{1}{W}\mathop{\sum }\limits_{j=1}^{n}\mathop{\sum }\limits_{r=1}^{W}{x}_{{k}_{r},j}{p}_{r}log{x}_{{k}_{r,}j}{p}_{r},$$\end{document} CDNFE i k = − 1 n 1 W ∑ j = 1 n ∑ r = 1 W x k r , j p r l o g x k r , j p r , with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${p}_{r}=\frac{{w}_{{in}}\left({g}^{k},{g}_{i,k}^{k}\right)}{{\sum }_{r=1}^{W}{w}_{{in}}\left({g}^{k},{g}_{i,k}^{k}\right)},$$\end{document} p r = w in g k , g i , k k ∑ r = 1 W w in g k , g i , k k , and 9 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{CDNFE}}_{o}^{k}=-\frac{1}{n}\frac{1}{L}\mathop{\sum }\limits_{j=1}^{n}\mathop{\sum }\limits_{z=1}^{L}{x}_{{k}_{z,}j}{p}_{r}log{x}_{{k}_{z,}j}{p}_{z,}$$\end{document} CDNFE o k = − 1 n 1 L ∑ j = 1 n ∑ z = 1 L x k z , j p r l o g x k z , j p z , with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${p}_{z}=\frac{{w}_{{out}}\left({g}^{k},{g}_{o,k}^{k}\right)}{{\sum }_{r=1}^{W}{w}_{{out}}\left({g}^{k},{g}_{o,k}^{k}\right)}$$\end{document} p z = w out g k , g o , k k ∑ r = 1 W w out g k , g o , k k \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${x}_{{k}_{r},j}$$\end{document} x k r , j and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${x}_{{k}_{z,}j}$$\end{document} x k z , j are the expression data of gene \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}_{r}^{k}$$\end{document} g r k and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}_{z}^{k}$$\end{document} g z k in the sample/cell \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j$$\end{document} j . [Step 6] Calculate the CDNFE score of the global causality directed network. At time point t, the CDNFE score for the perturbed sample is defined as follows, 10 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{CDNFE}}_{t}=\frac{{\sum }_{k=1}^{m}{{CDNFE}}_{t}^{k}}{m}$$\end{document} CDNFE t = ∑ k = 1 m CDNFE t k m The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{CDNFE}}_{t}$$\end{document} CDNFE t score at each time point t reflects the global perturbation caused by case samples/cells. An abrupt increase in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{CDNFE}}_{t}$$\end{document} CDNFE t score indicates that the point t may be a tipping point of critical transition. To evaluate the capacity of the CDNFE score in capturing critical state, we applied one-sample t -test ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S=\frac{\sqrt{n}\left({mean}(\hat{X})-x\right)}{{SD}\left(\hat{X}\right)}$$\end{document} S = n mean ( X ˆ ) − x SD X ˆ ), and the time point t can be conceived as critical state if the CDNFE score \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C{{DNFE}}_{t}$$\end{document} C DNFE t meets the following two conditions: (i) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{CD}{NFE}}_{t}$$\end{document} CD NFE t  >  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{CD}{NFE}}_{t-1}$$\end{document} CD NFE t − 1 , and (ii) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{CD}{NFE}}_{t}$$\end{document} CD NFE t is statistically different ( P  < 0.05) from the prior information. The top 5% genes in terms of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{CDNFE}}_{t}^{K}$$\end{document} CDNFE t K are DNBs in this work. The CDNFE method mainly characterizes the network fluctuations of molecules in a system rather than the differences or random fluctuations in the system, which is the core of quantifying the criticality or tipping point of the network, thus providing a robust and reliable early warning signal for the critical state. The detailed pseudocode is as follows: 1. Algorithm. Input: - Gene expression matrix - Parameters: λ, N Output: - Local CDNFE scores: CDNFE k t - Global CDNFE score: CDNFE t - Causal directed networks for each gene 2. Begin Algorithm: 3. 1. Co-expression Network Construction: 4. Apply WGCNA to GEM 5. Generate undirected co-expression network 6. 2. Directed Network Construction: 7. for each edge (gᵏ, gˡ) do: 8. if direction determination index W kl > 0 then 9. add directed edge: gᵏ → gˡ 10. end if 11. end for 12. for each gene gᵏ: 13. Calculate direction index: D k \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\leftarrow$$\end{document} ← λD k (in) + (1-λ)D k (out) 14. end for 15. for each directed edge (gᵏ → gˡ): 16. Calculate edge weight: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${w}_{t}\left({g}^{k},{g}_{j}^{k}\right)\leftarrow ln\left(\frac{\hat{\beta }}{\beta }\right)$$\end{document} w t g k , g j k ← l n β ˆ β 17. end for 18. 3. Local Network Extraction: 19. for each gene gᵏ: 20. Extract gᵏ-local in-degree network (edges pointing to gᵏ) 21. Extract gᵏ-local out-degree network (edges from gᵏ) 22. Combine to form Local Causal Directed Network for gᵏ 23. end for 24. 4. CDNFE Score Calculation: 25. for each gene gᵏ: 26. // In-degree network calculation 27. Compute CDNFEᵢᵏ for in-degree network 28. // Out-degree network calculation 29. Compute CDNFEₒᵏ for out-degree network 30. // Combine scores 31. W \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\leftarrow$$\end{document} ← number of in-degree neighbors, L \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\leftarrow$$\end{document} ← number of out-degree neighbors, N \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\leftarrow$$\end{document} ← W + L 32. CDNFEₜᵏ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\leftarrow$$\end{document} ← (W/N) CDNFE t k + (L/N) CDNFEₒᵏ 33. end for 34. 6. Global Score Calculation: 35. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{CDNFE}}_{t}\leftarrow \frac{{\sum }_{k=1}^{m}{{CDNFE}}_{t}^{k}}{m}$$\end{document} CDNFE t ← ∑ k = 1 m CDNFE t k m 36. Return out The sKLD method achieves early warning of critical disease states in single samples through the following steps: First, using a set of healthy samples as the reference background, it fits a Gaussian distribution to the expression values of each gene and constructs a reference probability distribution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P$$\end{document} P . Next, for an individual test sample, it calculates the cumulative area under the expression curves of each gene and constructs a perturbation distribution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q$$\end{document} Q . Finally, it computes the symmetric Kullback–Leibler divergence between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P$$\end{document} P and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q$$\end{document} Q to obtain the sKLD score. The core steps of the NIG method can be summarized as follows: First, using gene expression data from a set of reference samples and a single test sample, construct an individual-specific network through weighted gene co-expression network analysis; Next, extract gene module information from this network and construct local networks for each gene; Then, based on network flow entropy theory, calculate the average cumulative information gain for each central gene both when only reference samples are present and when the test sample is added; Finally, by comparing the difference between the two, obtain the NIG score quantifying the degree of disturbance to the global network caused by the individual sample.

To benchmark the sensitivity and robustness of our proposed method against established noise levels, we first applied CDNFE to a controlled numerical environment using a simulated gene regulatory network. A simulated dataset was generated using an 11-node regulatory network with positive or negative relationships (as depicted in Fig. 2A ) governed by a system of stochastic differential equations (File S1 ). This model was employed to identify the critical phase as the system approaches a bifurcation point. Gene regulatory network dynamics, including processes such as transcription and translation 10 , 11 , cyclic reactions 12 , nonlinear biological processes 11 , 13 , 14 , and other regulatory activities 11 , 13 , are commonly modeled using Michaelis–Menten kinetics or the Hill function 15 , 16 . By adjusting the parameter range from −0.50 to 0.25, a dataset was generated for numerical simulation. Fig. 2 Validation of CDNFE reliability and robustness using numerical simulations. A The 11-node network model used, with the DNB module highlighted in dashed lines. B Node-level score landscape showing coordinated entropy increases as the system approaches the critical state. C Schematic illustrating the CDNFE peak at the tipping point between stable states. D Global CDNFE score trajectory, accurately identifying the theoretical bifurcation point at control parameter p  = 0. E Performance benchmark against existing methods (NIG, KLD) under increasing noise levels. CDNFE provides a distinct, stable warning signal even at high noise intensities (right panel), outperforming other metrics. A The 11-node network model used, with the DNB module highlighted in dashed lines. B Node-level score landscape showing coordinated entropy increases as the system approaches the critical state. C Schematic illustrating the CDNFE peak at the tipping point between stable states. D Global CDNFE score trajectory, accurately identifying the theoretical bifurcation point at control parameter p  = 0. E Performance benchmark against existing methods (NIG, KLD) under increasing noise levels. CDNFE provides a distinct, stable warning signal even at high noise intensities (right panel), outperforming other metrics. To illustrate the contrasting dynamics more clearly between the normal and critical states, we presented the landscape of CNDFE scores for 11 nodes in Fig. 2B . As shown in Fig. 2 C, D, the CDNFE score for the 11-node network exhibited a significant change near the specific parameter value p  = 0, which indicates the tipping point or critical transition at a bifurcation point ( p  = 0). It was evident that when the system is far away from the tipping point, the CDNFE scores across all nodes remain stable and low. However, as the system nears the tipping point p  =  0 , some nodes exhibit markedly different behavior in terms of expression changes and network interactions, known as DNB, leading to a notable increase in the CDNFE score. This increase is indicative of the approaching tipping point or critical state. In addition, to demonstrate the robustness of the proposed method, we conducted a comparative analysis with the sample Kullback–Leibler divergence (sKLD) and Network Information Gain (NIG) methods 8 , 9 using data samples subjected to various levels of noise perturbation, as depicted in Fig. 2E . With the increase of noise strength, CDNFE method consistently delivers reliable early warning signals for the impending tipping point, with a more pronounced and sensitive score. In summary, these numerical simulations confirm that CDNFE is not only capable of detecting the theoretical bifurcation point with high precision but also demonstrates superior robustness compared to existing methods (sKLD and NIG). The method maintains reliable early-warning signals even under high-noise conditions, making it suitable for analyzing noisy biological datasets. Detailed simulation and calculation procedures are elaborated in the File S1 . Having validated the method in silico, we next sought to determine if CDNFE could identify a specific precancerous tipping point in real-world patient data derived from cervical tissue samples spanning multiple disease stages. The CDNFE method was implemented on cervical cancer data from Luohe Central Hospital ( n  = 9, spanning normal, precancerous, and cancer stages) to demonstrate its efficacy in pointing the critical state in tumors, particularly cervical cancer. For each individual Cervical squamous cell carcinoma (CESC) sample, we computed the individual CDNFE score. Figure 3A illustrates a significant rise in CDNFE scores as we transition from Normal to high-grade squamous intraepithelial lesions (HSIL) and from early CESC to late CESC ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P=\,0.0503$$\end{document} P = 0.0503 ). So, the pre-disease state of cervical cancer was identified at HSIL ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P=\,1.26E-06$$\end{document} P = 1.26 E − 06 ). The top 5% of genes with the highest CDNFE scores at the critical state were designated as DNBs. Figure 3B presents the landscape of local CDNFE scores, where a notable increase in CDNFE scores for DNBs is observed at HSIL. However, as shown in Fig. 3C , simply using the gene expressions of DNB genes cannot distinguish the critical state ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P=\,0.0680$$\end{document} P = 0.0680 at HSIL) from other states, but by using CDNFE scores of those DNB genes, we can detect the critical state as shown in Fig. 3A . To benchmark CDNFE against existing approaches, we applied NIG and sKLD to the Luohe dataset. Both methods detected a critical transition, albeit with notable differences. NIG identified the critical point at early CESC (Fig. S1a ), which deviates from the HSIL stage identified by CDNFE (Fig. S1b ). In contrast, sKLD aligned with CDNFE in identifying HSIL as the critical stage, although this alignment did not reach statistical significance ( P  = 0.6947). To assess the robustness of this finding, we expanded the discovery cohort by integrating our dataset with a public dataset ( GSE7803 ) after rigorous batch effect correction using ComBat (see “Methods”). When CDNFE was applied to the integrated cohort, the critical transition again occurred at the HSIL stage (Fig. S2 ). This consistency across independent cohorts confirms that the identified critical transition is not an artifact of small sample size and supports the generalizability of our findings. Fig. 3 The identification of the tipping point for cervical cancer based on CDNFE. A Global CDNFE scores across disease stages. The CDNFE scores are notably elevated during HSIL of CESC’s differentiation process compared to other stages, suggesting that HSIL represents a tipping point in the progression of CESC. B 3D landscape of the DNBs for CESC. The CDNFE score in DNBs increases significantly at HSIL, indicating the critical state. C Average expression of the same DNB genes, which fails to distinguish the critical state, demonstrating CDNFE’s superior sensitivity over standard expression analysis. D The dynamic evolution of the gene directed regulatory network constructed by signaling genes for this dataset. The nodes represent genes and directed edges represent the interactions between genes. The thickness of the edges indicates the magnitude of the causal weights. E A regulatory network centered on NECTIN4, with arrows indicating the causal direction between, and bar graphs showing the scores of genes at different stages. F , G Pseudotime trajectory inference. Trajectories based on raw gene expression ( F ) are dispersed, whereas those based on CDNFE scores ( G ) reveal a coherent developmental lineage with distinct bifurcation points. A Global CDNFE scores across disease stages. The CDNFE scores are notably elevated during HSIL of CESC’s differentiation process compared to other stages, suggesting that HSIL represents a tipping point in the progression of CESC. B 3D landscape of the DNBs for CESC. The CDNFE score in DNBs increases significantly at HSIL, indicating the critical state. C Average expression of the same DNB genes, which fails to distinguish the critical state, demonstrating CDNFE’s superior sensitivity over standard expression analysis. D The dynamic evolution of the gene directed regulatory network constructed by signaling genes for this dataset. The nodes represent genes and directed edges represent the interactions between genes. The thickness of the edges indicates the magnitude of the causal weights. E A regulatory network centered on NECTIN4, with arrows indicating the causal direction between, and bar graphs showing the scores of genes at different stages. F , G Pseudotime trajectory inference. Trajectories based on raw gene expression ( F ) are dispersed, whereas those based on CDNFE scores ( G ) reveal a coherent developmental lineage with distinct bifurcation points. The regulatory network, which integrates DNBs and their neighboring differentially expressed genes, is utilized to elucidate the network-level molecular regulatory mechanisms that drive the progression of cervical cancer. As depicted in Figs. 3D and S3 , the dynamic evolution of the Protein-Protein Interaction (PPI) directed network of CDNFE signal biomarkers. After the critical state, the gene expression patterns within the network changed significantly and the CDNFE scores showed a marked change from high to low or from low to high. The nodes symbolizing individual genes and the causal directed edges representing the interactions between these genes together form a dynamic causal directed network. This network structure helps to reveal the intricate regulatory mechanisms in the complex biological processes of cervical cancer and reveals the relationship between genes in many ways. Within these complex networks, we have identified some major sub-networks where the central genes are designated as DNBs. For example, in this key sub-network, NECTIN4 (as illustrated in Fig. 3E ) emerges as a central hub, with strong directed associations to other genes. Notably, NECTIN4 is highly expressed in cervical cancer, hinting at its potential role in disease progression 17 . The overexpression of NECTIN4 is known to trigger the PI3K/AKT signaling pathway 18 . ITGB6 is another gene implicated in this pathway. Figure 3E suggests that NECTIN4 may exert a regulatory influence on ITGB6 . Pseudotime analysis orders individual cells along a continuous trajectory based on similarities in molecular profiles, thereby approximating developmental or disease progression without requiring time-series data. Here, we applied trajectory inference to the same single-cell RNA-seq dataset using two alternative feature representations, one the conventional gene expression matrix (GEM), and the other a matrix of CDNFE scores computed for each gene in each cell. When pseudotime inference was performed using raw gene expression features (Fig. 3F ), cells exhibited a highly dispersed trajectory structure, characterized by broad branching and diffuse transitions between cell states. This dispersion reflects substantial transcriptional heterogeneity and stochastic noise at the single-gene expression level, which can obscure coordinated regulatory programs during critical state transitions. In contrast, trajectory inference based on CDNFE scores produced a more coherent and structured trajectory with three clear bifurcation points (Fig. 3G ). Unlike raw expression, CDNFE scores capture coordinated changes in network entropy among causally connected genes, thereby emphasizing collective regulatory instability rather than isolated expression fluctuations. As a result, CDNFE-based trajectories highlight synchronized shifts in cell state and reveal distinct lineage bifurcations associated with the precancerous tipping point. However, there were also commonalities. For instance, Smooth_muscle_cells and Keratinocytes were appearing in the same chronological order over time. Upon visualizing the key genes, we notably observed that the DNB gene FNDC3B exhibited high CDNFE scores across all cell types. It is known that FNDC3B is co-expressed with PRRX1 , ITGAV , and SKIL 19 , and in our findings, these genes are not only co-expressed but also appear to be regulated by FNDC3B . These neighboring genes also showed high CDNFE score in various cell types (Fig. S4 ). SKIL is recognized for facilitating cancer cells’ evasion of T-cell-mediated immunity via autophagy-driven suppression of the STING pathway 20 , and (Fig. S5 ) shows the expression of SKIL is low in T-cells (Fig. S6 ). Collectively, these results show that CDNFE identifies a tipping point at the HSIL stage characterized by a specific DNB module. This network-based signal successfully detects the transition that standard differential gene expression analysis fails to capture, highlighting the role of NECTIN4 and FNDC3B in the progression to malignancy. In many clinical practices and biological researches, differential expression genes are used to detect disease or mutant genes. However, we found that those genes without differential expression but sensitive to the CDNFE score have important functional roles in cancer progression, and are the candidates of “dark genes” 21 . To understand the biological mechanisms driving the identified tipping point, we next investigated whether the dark gene module uncovered by CDNFE is not only computationally detectable but also biologically and structurally meaningful. We analyzed the functional enrichment and spatial distribution of the DNB module and the “dark genes” within the tumor microenvironment (TME). We performed functional enrichment analyzes of DNBs to explore the functional correlation between DNB genes and cervical cancer progression (Table S1 ). The enrichment results showed that these DNBs were mainly enriched in the innate immune response, regulation of extrinsic apoptosis, cell adhesion molecule binding and other cancer-related signaling pathways (Fig. 4A ). Through literature searching, these pathways are indispensable in immune alterations, tumor cell proliferation, survival, and metastasis as the tumor evolves 22 , 23 . Moreover, to delve deeper into the potential biological functions of the “dark genes” and their adjacent first-order neighborhood differential genes, we conducted pathway enrichment analysis (Table S2 ). The process of identifying first-order neighborhood differential genes of dark genes is predicated on two fundamental criteria: their position as first-order neighbors of dark genes within the PPI network, and a statistically significant alteration ( P  < 0.05) in their gene expression levels before and after the critical state has been identified. As presented in the Fig. 4B , the pathways significantly enriched in “dark genes” were mainly focused on the virus life cycle, Chemical carcinogenesis and Interleukin-1 family signaling. It is widely acknowledged that infection with HPV is a significant etiological factor in cervical cancer development 24 . The integration of HPV, disruption of the E2 gene, usurpation of host PPIs, and enhancement of super-enhancers lead to the uncontrolled expression of oncoproteins, which promote cellular transformation and the progression to cervical neoplasia 25 . The activation of the viral life cycle pathway is likely to significantly accelerate HPV infection, ultimately contributing to the onset of cervical cancer. Similarly, under the influence of the chemical carcinogenesis pathway involving reactive oxygen species, mutations in purines and pyrimidines are induced, and chromatin proteins are oxidized. This triggers genomic instability, disrupts gene expression, and consequently promotes the initiation and progression of cancer 25 .IL1R2 plays a crucial role in the Interleukin-1 family signaling pathway by enhancing MHC-II expression on cancer-associated fibroblasts (CAFs), thereby increasing the number of regulatory T cells in the TME 26 . We can also extract from Figs. 4C and S9 that most of the 1st order differential genes are significantly enriched in multiple pathways. Significantly, the expression levels of first-order neighborhood differential genes exhibited shifts from high to low and vice versa, preceding and following the critical threshold state. The pathways in which first-order neighborhood differential genes are significantly enriched are also closely associated with cancer development. These pathways control cell proliferation and apoptosis 27 , cancer cell migration and invasion 28 , and immune modulation 29 . Fig. 4 The molecular signaling mechanism and spatial validation anchors the tipping point to tumor biology. A DNBs are involved in important biological processes and KEGG pathways for cancer. B ‘Dark genes’ are involved in important biological processes and KEGG pathways. The left side of the outer ring represents ‘dark gene’ detected by CDNFE, and the right side of the outer ring represents various biological processes in which these genes are involved. In the inner ring, the color and width of links indicate diverse enrichment pathways and significant levels of gene functions, respectively. C KEGG pathway and GO pathway enrichment analysis for the these1st-order DEG neighbors. The left side of this figure is a sankey representing the genes contained in each pathway, and the right side is a regular bubble plot, with the size of the bubbles indicating the number of genes to which the pathway belongs, and the color of the bubbles indicating the p -values. D Network controllability analysis shows a high proportion of dark genes are classified as driver nodes (orange and purple). E Dark genes have a significantly higher network out-degree than other genes, indicating greater regulatory influence (p = 0.034). F – H Spatial validation of the DNB signature in an independent invasive cervical cancer specimen. F DNB module score mapped to in-tissue Visium spots and overlaid on the matched histology image. G Tumor-enriched ROI overlay, where the ROI is defined from high epithelial and proliferation marker signal (top quantile of the composite marker score), shown on the same tissue section for spatial context. H DNB module scores are significantly higher in tumor-enriched ROI spots than in non-ROI spots (boxplot; Mann–Whitney U test). A DNBs are involved in important biological processes and KEGG pathways for cancer. B ‘Dark genes’ are involved in important biological processes and KEGG pathways. The left side of the outer ring represents ‘dark gene’ detected by CDNFE, and the right side of the outer ring represents various biological processes in which these genes are involved. In the inner ring, the color and width of links indicate diverse enrichment pathways and significant levels of gene functions, respectively. C KEGG pathway and GO pathway enrichment analysis for the these1st-order DEG neighbors. The left side of this figure is a sankey representing the genes contained in each pathway, and the right side is a regular bubble plot, with the size of the bubbles indicating the number of genes to which the pathway belongs, and the color of the bubbles indicating the p -values. D Network controllability analysis shows a high proportion of dark genes are classified as driver nodes (orange and purple). E Dark genes have a significantly higher network out-degree than other genes, indicating greater regulatory influence (p = 0.034). F – H Spatial validation of the DNB signature in an independent invasive cervical cancer specimen. F DNB module score mapped to in-tissue Visium spots and overlaid on the matched histology image. G Tumor-enriched ROI overlay, where the ROI is defined from high epithelial and proliferation marker signal (top quantile of the composite marker score), shown on the same tissue section for spatial context. H DNB module scores are significantly higher in tumor-enriched ROI spots than in non-ROI spots (boxplot; Mann–Whitney U test). Finally, we investigated whether the dark gene module uncovered by CDNFE is not only computationally detectable but also biologically and structurally meaningful. From a network perspective, driver node analysis revealed that dark genes disproportionately occupy critical structural positions within the regulatory network. Mathematically, out-degree represents the number of edges emanating from a node to others. Biologically, a gene with high out-degree functions as a ‘source’ node, propagating signals to numerous downstream targets. Many dark genes emerged as driver nodes in the minimal control set required to steer network dynamics (Fig. 4D ), and their out-degree was significantly higher compared to other genes ( p  = 0.034, Mann–Whitney U Test) 30 and are disproportionately positioned as “driver nodes” theoretically capable of controlling the network’s state (Fig. 4 D, E) 31 . To also provide definitive anatomical relevance, we sought to localize the DNB program within tissue using 10x Genomics Visium spatial transcriptomics 32 . We computed a DNB module score using the 75 DNB genes (63 detected in this dataset) and visualized the score across the section on the matched histology image (Fig. 4F ). We then define tumor-enriched regions based on tumor-associated signature (epithelial and proliferation marker signal) and tested whether DNB module scores were higher in tumor-enriched regions than in the remaining spots (Fig. 4G ). Quantitatively, DNB module scores were significantly enriched in tumor-enriched regions than outside (Fig. 4H ; Mann–Whitney U test). This directly links the abstract tipping point signal to the physical architecture of the tumor, proving that our computationally-derived signature marks the anatomical location of the disease. To address the limitation of a single cervical spatial section, we repeated the same scoring workflow in an independent Visium sample from head and neck squamous cell carcinoma ( GSM5494475 ) and observed the same qualitative enrichment of DNB signal in tumor-enriched regions (Fig. S13 ). A key feature of the tipping point was that it was not driven by conventional biomarkers. A volcano plot of the differential expression analysis between the tipping point and normal stages revealed that the vast majority of the DNB genes identified by CDNFE were not significantly changed—they were “dark” (Fig. 5A ). This was confirmed at the group level, where the average expression of the dark genes remained stable across all disease stages (Figs. 5B and S12 ), while the distribution of CDNFE scores shifted sharply, peaking precisely at HSIL (Fig. 5C ). Fig. 5 The tipping point is driven by a “dark gene” signature. A Volcano plot comparing HSIL vs. Normal. DNB genes (red dots) are predominantly located in the non-significant central region, defining them as “dark.” B The average expression of the DNB gene module is stable across disease progression. C In contrast, the distribution of CDNFE scores for the DNB genes shifts significantly, peaking at the HSIL tipping point. D Representative dark genes show stable expression (bar charts, left y -axis) but dynamic CDNFE scores (line plots, right y -axis) that peak at HSIL. E ROC curves for CDNFE scores of dark genes show strong predictive power for the tipping point. F ROC curves for the expression levels of the same genes perform no better than random chance. G The Area Under the Curve (AUC) values demonstrate that CDNFE scores vastly outper form expression levels as predictive metrics. H Regulation of related ‘dark genes’ and 1st order DEGs in the PI3K-AKT signaling pathway and MAPK signaling pathway. A Volcano plot comparing HSIL vs. Normal. DNB genes (red dots) are predominantly located in the non-significant central region, defining them as “dark.” B The average expression of the DNB gene module is stable across disease progression. C In contrast, the distribution of CDNFE scores for the DNB genes shifts significantly, peaking at the HSIL tipping point. D Representative dark genes show stable expression (bar charts, left y -axis) but dynamic CDNFE scores (line plots, right y -axis) that peak at HSIL. E ROC curves for CDNFE scores of dark genes show strong predictive power for the tipping point. F ROC curves for the expression levels of the same genes perform no better than random chance. G The Area Under the Curve (AUC) values demonstrate that CDNFE scores vastly outper form expression levels as predictive metrics. H Regulation of related ‘dark genes’ and 1st order DEGs in the PI3K-AKT signaling pathway and MAPK signaling pathway. Figure 5D illustrates a selection of “dark genes” that display minimal fluctuations in gene expression levels but exhibit substantial changes in CDNFE values. Additional “dark genes” relevant to this dataset are detailed in Table S3 . Moreover, these “dark genes” have significant functional roles in cancer processes. For example, DNAJB4 has been implicated in the development and progression of various cancers. Regarding cancer, a growing body of research indicates that DNAJB4 acts as a tumor suppressor, which can inhibit cancer cell proliferation, invasion, and metastasis. The high expression of DNAJB4 in cancer cells can suppress their proliferation and tumorigenic potential, as well as reduce cell mortality and invasiveness. This is achieved by downregulating cyclin D1 , upregulating STAT1 and p21WAF1 , and subsequently activating the STAT1 signaling pathway 33 . Meanwhile, DNAJB4 has been identified as a potential promoter of metastasis 34 . ASPM plays an important role in the formation of the mitotic spindle during cell division, initiating the assembly of microtubules at the centrosomes. Extensive research has established a connection between ASPM and tumor progression, as well as its association with poor survival outcomes 35 . Additionally, LINC01133 has been shown to enhance the proliferation, migration, and epithelial-mesenchymal transition (EMT) in cervical cancer (CC) cells. LINC01133 competes with miR-4784 for binding to AHDC1 , leading to increased AHDC1 expression and intensifying the malignant phenotypes and EMT process in CC cells 36 . To prove the utility of this concept, we tested whether CDNFE scores or raw expression values were better predictors of the tipping point. Using Receiver Operating Characteristic (ROC) analysis 37 , we demonstrated that the CDNFE scores of dark genes had vastly superior predictive power, with AUCs often exceeding 0.7 (Fig. 5 E, F). By contrast, the AUCs for their near-random expression levels were often below 0.5 (Fig. 5G ). This demonstrates that the CDNFE score captures hidden regulatory signals that define the state transition, establishing dark genes as a novel class of biomarkers beyond the reach of differential expression. We also performed the functional analysis for ‘dark genes’ and their first-order neighborhood differential genes, which shows that some genes are enriched into the key biological processes or functions associated with the development of cervical cancer. We focused on the analysis of the KEGG pathway on the PI3K-AKT signaling pathway and MAPK signaling pathway that were most relevant to the progression of cervical cancer. The Rac1 activated by the “dark gene” NECTIN4 was involved in the PI3K-AKT signaling pathway (Fig. 5H ). The 1st order DEGs, such as ITGBV , ITGB6 and GNB1 , activated the PI3K-AKT signaling pathway. After the critical state (HSIL), signaling genes ( NECTIN4 ) were highly expressed, which might have triggered the phosphorylation of PI3K and AKT proteins and further reduced the expression of the apoptosis inhibitor GSK3B , and finally enhanced the cell cycle. PI3K-AKT signaling pathway plays a crucial role in tumor growth, cancer metabolism and metastasis 38 . Proteins encoded by the “dark genes” MAP3K20 and 1st order DEGs ( CACNA2D3 , PPP3CA , MAPK10 , HSPA6 ) are involved in the MAPK signaling pathway, which is a potential pathway for the development of cervical cancer (Fig. 5H ). MKK7 activity can be enhanced through either MKK7 autophosphorylation or phosphorylation by upstream MAP3K20 in the kinase domain. Upon reaching a critical state, the expression of MAP3K20 was notably upregulated. Concurrently, the expression of MAPK10 , regulated by MKK7 , also saw a significant increase. Additionally, HSPA6 has a negative regulatory effect on MAPK10 , high expression of HSPA6 inhibits the increase in MAPK10 levels. Furthermore, the elevation of intracellular calcium ( Ca 2+ ) induced by CACNA2D3 indirectly impacts PPP3CA , which in turn affects NFTC1 . The upregulation of NFTC1 proteins enhances tumor invasiveness and promotes tumor progression 39 . As the same time, tumor growth depends upon an adequate supply of oxygen and nutrients, which are delivered via blood vessels. Angiogenesis is a critical mediator of tumor cell metastasis. Within the TME, the JAM and EPHA pathways (Fig. S8 ) are instrumental in modulating MAPK signaling by governing the interactions between MAPK and other proteins, which affects cell proliferation and survival 40 , 41 . In the EPHA pathways, EFNA1 and its receptor EPHA2 (Fig. S9 ) are implicated in angiogenesis, interacting with other molecules such as Vav and PI3K, thereby facilitating the progression of cervical cancer 42 – 44 . These multi-modal validations provide a biological anchor for our computational findings. We demonstrate that the “dark gene” signature is not a statistical artifact but a biologically coherent module enriched in cancer-driving pathways (PI3K-AKT, MAPK) and structurally positioned as network hubs. Crucially, the spatial localization of this signature to malignant epithelial nests confirms its direct relevance to the physical architecture of the tumor. Finally, to determine whether these identified dark genes are merely passengers or essential drivers of cancer cell viability, we integrated proteomic data and large-scale CRISPR-Cas9 dependency screens to assess their functional necessity. Proteogenomic analysis was utilized to confirm that “dark genes” are indeed translated into proteins and may exert functional influences in the progression of cervical cancer, even in the absence of differential expression patterns. Supporting proteomic evidence pertaining to these “dark genes” was derived from pertinent research studies 45 . For all 30 dark genes of CESC, the result shows that 27 proteins corresponding to the “dark genes” are identified with high confidence by isobaric tandem mass tags-based proteomic analysis, in which six differentially expressed proteins corresponding to the “dark genes” are identified with p -value ≤ 0.05 and log2.median.fold-change ≥1 or ≤ −1 46 . Especially, NCK Adaptor Protein 1, corresponding to the dark gene NCK1 was found to may play important roles in cancer angiogenesis 47 . In addition, in the tumors, six proteins, that are Glutathione S-Transferase Mu 3, Assembly Factor For Spindle Microtubules, DnaJ Heat Shock Protein Family (Hsp40) Member B4, OCIA Domain Containing 2, Immunoglobulin Kappa Variable 3-20 and Cornulin, corresponding to the dark genes GSTM3 , ASPM , DNAJB4 , OCIAD2 , IGKV3-20 and CRNN , are differentially expressed (≥20-fold) relative to the paired NATs. Assembly Factor For Spindle Microtubules and Glutathione S-Transferase Mu 3 are upregulated and have a strong association with survival of patients 48 . On the other hand, by significance test with hypergeo-metric distribution for the results of CESC, the probability of 27 proteins with high confidence among any 30 genes (randomly chosen from 1000 genes) is 0.1097586, and further the probability of six differential proteins with p -value ≤ 0.05 and log2.median. fold-change ≥1 or ≤ −1 expression from any 30 genes (randomly chosen from 1000 genes) is 0.1799619, which is clearly far significant. Therefore, our results of the 27 proteins with high confidence and the 6 differential expression proteins of our 34 dark genes at the protein level are much significant, comparing with other genes, which showed the functional roles of the dark genes at the protein level. Having established a predictive dark gene signature, we initiated our multi-modal validation framework to confirm its biological significance. First, we tested if these genes were functionally essential for the resulting cancer state (Fig. 6A ). By interrogating the DepMap CRISPR screen database across 18 cervical cancer cell lines 49 , we found that a remarkable 73.3% (22 of 30) of testable dark genes were essential for survival in at least one cell line (Fig. 6 B, C). To illustrate this point, we examined SARS1 , a representative dark gene. Although it is not differentially expressed and thus overlooked by conventional methods, SARS1 showed strong essentiality across diverse cancer types (Fig. 6D ), underscoring those dark genes can encode core fitness functions hidden from static expression-based analysis. Fig. 6 Functional essentiality validation of the dark gene signature. A Validation strategy: Dark genes predicted by CDNFE are queried in the DepMap CRISPR database to confirm functional essentiality in cancer cells. B Dependency Landscape: Heatmap shows the dependency scores for dark genes across 18 cervical cancer cell lines. Negative scores (blue) indicate that gene knockout leads to reduced cell viability (essentiality). C Essentiality Rate: Donut chart summarizing that 73.3% of testable dark genes are essential for survival in cervical cancer cells, confirming they play critical biological roles despite lacking differential expression. D Pan-Cancer Validation: Dependency profile of the top hit, SARS1, is shown to be a pan-cancer dependency, indicating its fundamental role in cell survival. A Validation strategy: Dark genes predicted by CDNFE are queried in the DepMap CRISPR database to confirm functional essentiality in cancer cells. B Dependency Landscape: Heatmap shows the dependency scores for dark genes across 18 cervical cancer cell lines. Negative scores (blue) indicate that gene knockout leads to reduced cell viability (essentiality). C Essentiality Rate: Donut chart summarizing that 73.3% of testable dark genes are essential for survival in cervical cancer cells, confirming they play critical biological roles despite lacking differential expression. D Pan-Cancer Validation: Dependency profile of the top hit, SARS1, is shown to be a pan-cancer dependency, indicating its fundamental role in cell survival. In conclusion, this functional interrogation confirms that the dark genes identified by CDNFE are fundamental to cancer cell survival. The high rate of essentiality (73.3%) observed in CRISPR screens, combined with proteomic validation, strongly suggests that CDNFE successfully uncovers hidden therapeutic targets that would be discarded by traditional differential expression analyzes.

Detecting early warning signals of critical transitions in disease progression remains a major challenge in systems biology and translational medicine. Phenotypic similarities and modest changes in mean expression between pre-transition and critical states often obscure the onset of disease tipping points, limiting opportunities for early intervention. In this research, we introduce CDNFE, a model-free, data-driven method that leverages directed network dynamics to identify critical states. Unlike conventional biomarker strategies that rely on differential expression, CDNFE captures subtle but coordinated network-level changes, enabling the detection of hidden regulators that conventional analyzes overlook. Our study establishes the clinical context by accurately identifying the HPV-driven progression from normal tissue to the precancerous tipping point (HSIL). Beyond detection, our findings elucidate a detailed molecular mechanism driving this transition. We show that stromal signals (TGF-β, ECM) from CAFs are integrated by the key discovery genes, FNDC3B and NECTIN4. This integration appears to trigger the PI3K/AKT signaling pathway, promoting EMT and invasion within a broader context of immune evasion and angiogenesis (Fig. 7B ). This aligns with existing literature suggesting that NECTIN4 overexpression activates PI3K/AKT signaling to enhance cell proliferation and migration 18 , while FNDC3B has been implicated in mediating EMT in other carcinomas 15 . The ability of CDNFE to position these genes upstream of such critical pathways highlights its utility in decoding the regulatory logic of disease progression. Fig. 7 CDNFE-based identification of tipping-point biology in cervical cancer. A Schematic overview of histological progression from normal epithelium through precancerous lesions to invasive cervical cancer, overlaid with characteristic CDNFE score dynamics. CDNFE scores remain low in normal tissue, spike sharply at the precancerous tipping point, and stabilize at later cancer stages. B Proposed microenvironmental signaling landscape at the tipping point, integrating CDNFE-identified genes and pathways. Numbered interactions indicate an inferred sequence of key biological events: (1) Cancer-associated fibroblasts (CAFs) secrete TGF-β into the local microenvironment. (2) TGF-β signaling promotes HSIL progression and epithelial remodeling. (3) Extracellular matrix (ECM) remodeling facilitates tumor–stroma interactions and epithelial invasion. (4) CDNFE-identified “dark genes” (e.g., FNDC3B, NECTIN4) activate PI3K–AKT signaling. (5) PD-L1 expression on epithelial cells engages PD-1 on T cells, contributing to immune suppression. (6) Downstream activation of angiogenesis (VEGF) and epithelial–mesenchymal transition (EMT) promotes malignant progression. C The clinical relevance of the dark gene signature is supported through three validation pillars: Structural Validation identifies dark genes as key driver nodes within the regulatory network; Functional Validation using DepMap CRISPR screens reveals that 73.3% of these genes are essential for cancer cell survival; and Anatomical Validation using spatial transcriptomics confirms their specific enrichment within malignant tumor nests. A Schematic overview of histological progression from normal epithelium through precancerous lesions to invasive cervical cancer, overlaid with characteristic CDNFE score dynamics. CDNFE scores remain low in normal tissue, spike sharply at the precancerous tipping point, and stabilize at later cancer stages. B Proposed microenvironmental signaling landscape at the tipping point, integrating CDNFE-identified genes and pathways. Numbered interactions indicate an inferred sequence of key biological events: (1) Cancer-associated fibroblasts (CAFs) secrete TGF-β into the local microenvironment. (2) TGF-β signaling promotes HSIL progression and epithelial remodeling. (3) Extracellular matrix (ECM) remodeling facilitates tumor–stroma interactions and epithelial invasion. (4) CDNFE-identified “dark genes” (e.g., FNDC3B, NECTIN4) activate PI3K–AKT signaling. (5) PD-L1 expression on epithelial cells engages PD-1 on T cells, contributing to immune suppression. (6) Downstream activation of angiogenesis (VEGF) and epithelial–mesenchymal transition (EMT) promotes malignant progression. C The clinical relevance of the dark gene signature is supported through three validation pillars: Structural Validation identifies dark genes as key driver nodes within the regulatory network; Functional Validation using DepMap CRISPR screens reveals that 73.3% of these genes are essential for cancer cell survival; and Anatomical Validation using spatial transcriptomics confirms their specific enrichment within malignant tumor nests. CDNFE differs from traditional approaches by integrating causality-directed network structure with entropy dynamics, rather than relying solely on static differential expression. This allows it to resolve the complexity of biological interactions and highlight markers that signal instability at the systems level. Applied across simulations and two independent cervical cancer datasets, CDNFE consistently pinpointed the critical transition, demonstrating both robustness and generalizability. Applied to cervical cancer, CDNFE consistently identified a tipping point at the precancerous (HSIL) stage, both in our single-cell discovery cohort and in an independent bulk transcriptome dataset. A central contribution of this work is the conceptualization and validation of “dark genes” - biomarkers whose expression levels remain stable but whose entropy-based scores sharply peak at the critical state. This property makes them invisible to traditional DEG-based analyzes but central to the transition itself. Structurally, we found that dark genes are disproportionately represented as driver nodes in regulatory networks, suggesting they exert control over system dynamics through network rewiring rather than abundance changes. This offers a theoretical explanation for why some therapeutic targets fail to show overexpression; their pathogenicity lies in their connectivity, not their quantity. To translate these computational findings into biological relevance, we employed a multi-modal validation framework. Spatially, we bridged the gap between abstract network scores and tumor architecture using spatial transcriptomics. The localization of the DNB dark gene signature to tumor-enriched regions confirms that these signals are intrinsic to the tumor cells driving the disease. Functionally, the DepMap analysis provided critical evidence of clinical utility. The finding that 73.3% of the identified dark genes are essential for cervical cancer cell survival suggests that these network-derived markers are not merely passengers but are fundamental to tumor viability. Finally, pathway enrichment analysis revealed associations with core cancer processes, including PI3K-AKT signaling pathway. This positions dark genes, particularly FNDC3B and NECTIN4 , as actionable therapeutic targets for precision oncology. Together, these results suggest that FNDC3B and NECTIN4 may function as drivers of precancer progression through PI3K/AKT‑mediated EMT. To further contextualize the performance and biological output of CDNFE, we additionally compared the identified gene sets with those derived from DNFE 50 and TNFE 51 on the same dataset ( GSE63514 ) (Fig. S14 ). A subset of genes, including CLCA2 and RPL11, was consistently identified across all three methods. In contrast, CDNFE uniquely prioritized genes with established roles in EMT and cell migration, such as MSMO1, CDKN1B, HADHB, RPL14, and DEFB1, which were either absent or substantially lower ranked by the other methods. These findings suggest that by integrating causal directionality and second‑order neighborhood analysis, CDNFE may offer an advantage in uncovering regulatory relationships with stronger mechanistic relevance—identifying genes that appear to exert notable influence on network dynamics and potentially play important roles in disease progression. Despite these promising results, our study has limitations. The discovery cohort was modest in size ( N  = 9), and although findings were validated across independent datasets ( GSE63514 ), larger prospective cohorts will be required to fully characterize the heterogeneity of the tipping point across diverse patient populations. Moreover, while CRISPR screens and network analysis provide functional support, experimental validation in cervical cancer models will be critical to confirm causality. In conclusion, we demonstrate that CDNFE provides a robust framework for detecting tipping points and reveals dark genes as a novel class of network biomarkers. By integrating entropy dynamics with causal network structure, CDNFE advances our ability to pinpoint hidden regulators of disease progression. Future work extending CDNFE to additional cancers and integrating with single-cell multi-omics and spatial data holds promise for establishing dark genes as clinically relevant biomarkers for early detection and therapeutic targeting in cervical cancer and beyond.

Cervical cancer ranks as the fourth most prevalent cancer among women, posing a significant threat to global women’s health. Global clinical data projects that by 2020, there will be approximately 600,000 new cases of cervical cancer, with an estimated 340,000 fatalities 1 . The early signs of cervical cancer are often subtle or absent, which complicates the detection process, allowing the tumor to advance and infiltrate deeper into the tissues before it is diagnosed 2 . The World Health Organization has consistently emphasized that cervical cancer is a preventable disease. When detected at an early stage and treated promptly, it can be effectively cured 3 . Therefore, pinpointing the critical tipping point and identifying the key biomarkers associated with this transition is crucial for enhancing the cure rate of cervical cancer. The precancerous phase marks the boundary between a relatively healthy state and a malignant one, underscoring its vital role in the progression of the disease. When cervical cancer passes through the precancerous state, it can rapidly progress to the cancerous state. Unlike the irreversible state of cancer, identifying the precancerous state offers the potential to prevent the further spread of cancer cells and to control the progression of the disease through timely interventions 4 . Significant efforts have been invested in the discovery of biomarkers to improve the diagnosis of the post-transitional phase. Traditional approaches based on differential gene expression have identified markers such as MCM2 , TOP2A , and BLM 4 . However, it remains difficult to detect the critical state because gene expression may not differ much between the normal state and the precancerous state. To identify critical signals that precede transitions in biological systems, the theory of dynamic network biomarkers (DNBs) 5 has been put forth. DNBs are composed of a group of molecules, genes, or proteins which has strong correlation with each other. Unlike the traditional biomarkers, the concentrations of DNB molecules are collectively fluctuated near the critical state, rather than constant values or random fluctuations of molecules. Based on the DNB theory, a multitude of techniques have been devised and effectively utilized in the detection of critical conditions associated with intricate diseases and biological phenomena 6 , 7 . However, in real applications, it is difficult to obtain multiple samples for everyone, which limits the utilization of DNB theory in biological research and clinical practice. Moreover, the majority of these calculational approaches focus on identifying correlations between molecules or genes through undirected networks, frequently overlooking the valuable information embedded within the underlying cellular populations. Conversely, causal directed networks can be constructed to infer the directionality and potential causal relationships among genes or proteins. This is essential for hypothesizing regulatory mechanisms, gaining a deeper comprehension of potential regulatory interplay, and pinpointing candidate key regulatory elements and pathways. Therefore, it is challenging to develop an effective and robust method to detect the critical state of cancer based on an individual basis. Here, we proposed a novel computational method, causality directed network flow entropy (CDNFE) (Fig. 1 ), to predict early warning signals before complex diseases exacerbate. In contrast, our approach concentrates on the second-order neighbors that have causal direct interactions with any of the first-order neighbors. Firstly, we constructed a specific directed network at a particular time point by utilizing the WGCNA network and employing a direction determination index. Then, we built causal directed networks based on causal intensity metrics. Finally, we calculated the local CDNFE score for each local causality-directed network, leveraging the information from the causality-directed network. A significant rise in CDNFE is indicative of an approaching critical state or a pre-deterioration condition (Fig. 1 ). Our proposed CDNFE method is an effective tool for detecting critical states of complex biological processes, which has the following advantages: (i) We applied our methods to a real dataset to determine the critical state of cervical cancer. The samples were extracted from the patients’ cervical tissues, encompassing a spectrum of conditions from human papillomavirus (HPV) infection to precancerous lesions and cervical cancer. (ii) In the CDNFE method, we concentrated on second-order neighbors, which more accurately capture the structure of networks, thereby enhancing validity and simplifying the network complexity. (iii) The CDNFE method was based on causality-directed networks. This innovative foundation enabled us to gain a deeper comprehension of the regulatory interactions among genes and delve into the analysis of gene regulatory networks. We discovered that in the gene regulatory network, in the PI3K/AKT signaling pathway, FNDC3B and NECTIN4 as representative dark gene regulators of cervical cancer progression. (iv) We demonstrate that the CDNFE-derived dark gene signature is (1) essential for cancer cell survival, (2) structurally critical as driver nodes within the regulatory network, and (3) spatially localized to malignant epithelial nests. Fig. 1 Schematic diagram of the CDNFE method. a Data Collection: Samples representing different stages of cervical cancer progression are collected for scRNA sequencing. b Network Construction: A gene expression matrix is used to build a co-expression network, which is then converted into a directed network using an index that infers regulatory directionality. This network is refined into local causality-directed networks. c CDNFE Score Calculation: For each local network, an entropy-based score is calculated for incoming (in-degree) and outgoing (out-degree) information flow. These are combined to generate a local CDNFE score, and the average of all local scores gives the global CDNFE score for the sample. a Data Collection: Samples representing different stages of cervical cancer progression are collected for scRNA sequencing. b Network Construction: A gene expression matrix is used to build a co-expression network, which is then converted into a directed network using an index that infers regulatory directionality. This network is refined into local causality-directed networks. c CDNFE Score Calculation: For each local network, an entropy-based score is calculated for incoming (in-degree) and outgoing (out-degree) information flow. These are combined to generate a local CDNFE score, and the average of all local scores gives the global CDNFE score for the sample. In order to verify the robustness and validity of the CDNFE method, we performed numerical simulations based on datasets generated by the artificial gene regulatory network. These simulations demonstrated consistent stability and robustness in capturing tipping points. Notably, as the strength of noise increased, CDNFE still performed better than existing methods 8 , 9 in detecting early-warning signals and identifying DNBs. The CDNFE method also was applied to a public cervical dataset ( GSE63514 ) (File S4 and Figs. S7 and 8 ) and the critical state identified were consistent with the clinical samples from Luohe Central Hospital. The key biomarkers were also verified by a series of functional enrichment analyzes.

Supplementary Information Supplementary Information