BCTI: a Bayesian network-based method for revealing critical transitions in complex biological systems

Bayesian network fundamentals

A Bayesian network is a probabilistic graphical model that represents conditional dependencies among variables (e.g., gene expression levels) using a directed acyclic graph (DAG). In this framework, nodes correspond to genes, and directed edges represent potential regulatory relationships. The goal of Bayesian network structure learning is to identify the optimal network $\tilde{G}$ that best explains the observed data $D$, i.e., the network that maximizes the posterior probability $P\left(G|D\right)$ (Fig. 1B). According to Bayes’ theorem, the posterior probability of a network $G$ given data $D$ is:

(1) $$\begin{array}{c}P\left(G|D\right)\propto P\left(D|G\right)P\left(G\right)\end{array}.$$Notably, maximizing the posterior probability is equivalent to maximizing the product of the likelihood $P\left(D|G\right)$ and the prior $P\left(G\right)$:

(2) $$\begin{array}{c}\tilde{G}=\underset{G}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}P\left(G|D\right)=\underset{G}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}P\left(D|G\right)P\left(G\right).\end{array}$$

For a given Bayesian network structure $G$, the likelihood can be factorized as:

(3) $$\begin{array}{c}P\left(D|G\right)=\prod _{s=1}^{m}P\left(X_s|Pa\left(X_s\right)\right)=\prod _{s=1}^m\prod _{j=1}^{n}P\left(x_{sj}|Pa\left(X_s\right)\right)\end{array}$$where $X_s$ represents the random variable corresponding to the expression level of gene $g_s$, and $Pa\left(X_s\right)$ represents the parent nodes (regulators) of $X_s$ in $G$. The term $x_{sj}$ refers to the observed expression level of gene $g_s$ in sample $j$. This factorization leverages the conditional independence properties encoded in the DAG, allowing efficient modeling of high-dimensional data.

Bayesian networks provide a powerful framework for modeling the causal relationships among genes. However, traditional Bayesian network models are primarily designed for discrete variables, whereas gene expression data are inherently continuous. To address this mismatch, we adopt structural equation models (SEMs) as continuous-variable counterparts within the Bayesian framework.

SEMs model each gene as a linear combination of its parent genes and a random noise term. This formulation allows us to leverage the causal interpretability of Bayesian networks while accommodating the continuous nature of high-throughput transcriptomic data. Moreover, SEMs allow for the estimation of model fit via residuals, which directly feed into our network scoring criterion.

For continuous gene expression data, the conditional probabilities $P\left(X_s|Pa\left(X_s\right)\right)$ can be modeled using SEMs, which represent each variable $X_s$ as a linear combination of its parent nodes $Pa\left(X_s\right)$ with an added noise term:

(4) $$\begin{array}{c}{x}_s={\mathit{\theta }}_s{x}_{Pa\left(X_s\right)}+{ϵ}_s,\end{array}$$where ${\mathit{x}}_s\in {\mathbb{R}}^{1\times n}$ is the expression vector of $X_s$, ${\mathit{x}}_{Pa\left(X_s\right)}\in {\mathbb{R}}^{p\times n}$ represents the expression matrix of $Pa\left(X_s\right)$, where $p$ is the number of the parents, ${\mathit{\theta }}_s\in {\mathbb{R}}^{1\times p}$ is the regression coefficient vector to be estimated, and ${\mathit{ϵ}}_s\in {\mathbb{R}}^{1\times n}$ is a Gaussian noise vector.

Furthermore, we develop a network score, i.e., $H$ score, based on the minimum description length (MDL) principle in information theory to guide the Bayesian network structure learning, which considers both the goodness-of-fit to the data and the network’s complexity:

(5) $$\begin{array}{c}H\left(G\right)=\sum _{s=1}^m\left[F\left(X_s,Pa\left(X_s\right),{\hat{\mathit{\theta }}}_s^{mle}\right)+\lambda \left|{\hat{\mathit{\theta }}}_s^{mle}\right|\log n\right],\end{array}$$where:

(6) $$\begin{array}{c}F\left(X_s,Pa\left(X_s\right),{\hat{\mathit{\theta }}}_s^{mle}\right)=\sum _{v=1}^n{\left[x_{sv}-{\hat{\mathit{\theta }}}_s^{mle}{\mathit{x}}_{Pa\left(X_s\right)v}\right]}^2,\end{array}$$

${\mathit{x}}_{Pa\left(X_s\right)v}\in {\mathbb{R}}^{p\times 1}$ represents the expression vector of the parents $Pa\left(X_s\right)$ in the $v$th sample, ${\hat{\mathit{\theta }}}_s^{mle}=\underset{{\theta }_s}{\arg min}F\left(X_s,Pa\left(X_s\right),{\mathit{\theta }}_s\right)\in {\mathbb{R}}^{1\times p}$ is the maximum likelihood estimate of the parameter ${\theta }_s$ using least squares when $F\left(X_s,Pa\left(X_s\right),{\mathit{\theta }}_s\right)$ is the minimum, and $\left|{\hat{\mathit{\theta }}}_s^{mle}\right|$ is the number of to-be-estimated parameters in ${\hat{\mathit{\theta }}}_s^{mle}$ (in this study, it is the number of parents, i.e., $\left|{\hat{\mathit{\theta }}}_s^{mle}\right|=\left|\mathit{P}\mathit{a}\left(X_s\right)\right|$). Noteworthily, $\left[F\left(X_s,Pa\left(X_s\right),{\hat{\theta }}_s^{mle}\right)+\lambda \left|{\hat{\theta }}_s^{mle}\right|\log n\right]$ represents the evaluation score for the local subnetwork consisting of $X_s$ and its parents, called a local subnetwork score, while $H\left(G\right)$ represents the sum of the local subnetwork scores, with each subnetwork corresponding to a specific gene in $G$. Here, $F\left(X_s,\mathit{P}\mathit{a}\left(X_s\right),{\hat{\mathit{\theta }}}_s^{mle}\right)$ measures the goodness of fit to the data, and $\lambda \left|{\hat{\mathit{\theta }}}_s^{mle}\right|\log n$ quantifies the network complexity. The regularization parameter $\lambda$ controls the trade-off between the goodness of fit and the network complexity.

Workflow of BCTI

Given gene expression data with $m$ genes from different time points, to investigate the dynamic changes in gene regulatory networks and uncover critical states in complex biological systems, the algorithm of BCTI is proposed to simultaneously perform quantitative and qualitative analyses of system states from a network-level perspective:

Step 1. Construct the undirected correlation basis network (CBN) at each time point.

In this study, mutual information is a measure used to quantify the correlation between two random variables, capturing both linear and nonlinear dependencies (Kraskov, Stögbauer & Grassberger, 2004). For random variables $X_i$ and $X_j$, $i,j\in \left\{1,2,\dots ,m\right\}$, the mutual information between $X_i$ and $X_j$ is defined as the following form:

(7) $$\begin{array}{c}I\left(X_i,X_j\right)=\sum _{x\in X_i,y\in X_j}p\left(x,y\right)\mathrm{l}\mathrm{o}\mathrm{g}{\displaystyle \frac{p\left(x,y\right)}{p\left(x\right)p\left(y\right)}},\end{array}$$where $p\left(x\right)$ (or $p\left(y\right)$) represents the probability that the value of the variable $X_i$ (or $X_j$) is equal to $x$ (or $y$), and $p\left(x,y\right)$ represents the joint probability of $X_i$ and $X_j$. To improve computational stability and efficiency, we followed a previous study (Zhao et al., 2016) and simplified Eq. (7) under the assumption that the variables follow Gaussian distributions:

(8) $$\begin{array}{c}I\left(X_i,X_j\right)=-{\displaystyle \frac{1}{2}\log \left(1-{\rho }_{ij}^2\right)},\end{array}$$where ${\rho }_{ij}$ is the Pearson correlation coefficient between $X_i$ and $X_j$. To prevent numerical instability when ${\rho }_{ij}$ approaches 1, a small constant ( $\epsilon ={10}^{-5}$) was added inside the logarithmic term during the computation:

(9) $$\begin{array}{c}I\left(X_i,X_j\right)=-{\displaystyle \frac{1}{2}\log \left(1-{\rho }_{ij}^2+\epsilon \right).}\end{array}$$

More computational details are provided in Note S2. To construct a CBN that captures adjacency relationships among genes, we set a default threshold $T=0.1$ to evaluate the relevance between genes, following previous network-based studies (Chaitankar et al., 2009, 2010; Yang et al., 2016; Zhang et al., 2012), where a similar cutoff was used to retain potentially relevant associations while excluding unlikely connections. Specifically, if $I\left(X_i,X_j\right)\ge T$, $X_i$ is deemed related to $X_j$; otherwise, if $I\left(X_i,X_j\right)<T$, the mutual information is insufficient to establish relevance, thus $X_i$ and $X_j$ are treated as independent (Liu et al., 2016), $i,j\in \left\{1,2,\dots ,m\right\}$. Essentially, the CBN illustrates the potential regulatory relationships between genes but does not indicate the direction of regulation.

Step 2. Randomly generate an initial DAG.

To enhance the likelihood of identifying the global optimal GRN, we repeatedly generate different initial DAGs from the undirected CBN structure. Taking the $k$th repetition as an example ( $k\in \left\{1,2,\dots ,I\right\}$), an initial DAG, conforming to the CBN structure, is randomly generated and denoted as $G_{k0}$.

Step 3. Determine a potential optimal network ( ${\tilde{G}}_k$) using a hill-climbing search.

Letting ${\tilde{G}}_k=G_{k0}$, we first apply the add-edge, delete-edge, and reverse-edge operations to $G_{k0}$, generating three sets of adjusted networks, denoted as $C_1$, $C_2$, and $C_3$, respectively (Note: these operations only target the edges identified in the CBN). Next, the network with the lowest $H$ score in each of the sets, called an alternative network, is identified for subsequent analysis and is denoted as follows:

(10) $$\begin{array}{c}{G}_{k1}=\underset{G_{k1}^i\in C_1}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\left\{H\left(G_{k1}^i\right)\right\},\end{array}$$

(11) $$\begin{array}{c}{G}_{k2}=\underset{G_{k2}^j\in C_2}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\left\{H\left(G_{k2}^j\right)\right\},\end{array}$$

(12) $$\begin{array}{c}{G}_{k3}=\underset{G_{k3}^q\in C_3}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\left\{H\left(G_{k3}^q\right)\right\}.\end{array}$$

Subsequently, the $H$ scores of these alternative networks are compared to obtain a temporary candidate network, denoted as $G_k$:

(13) $$\begin{array}{c}{G}_k=\underset{G_{kt}\in \left\{G_{k1},G_{k2},G_{k3}\right\}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\{H(G_{kt})\}.\end{array}$$If the $H$ score of $G_k$ is lower than that of ${\tilde{G}}_k$, then ${\tilde{G}}_k=G_k$. Next, we set ${\tilde{G}}_k$ as the new initial directed network and continue performing edge adjustments, applying the adjustment that yields the minimum $H$ score to refine ${\tilde{G}}_k$ and produce a new $G_k$. This process is repeated iteratively until the $H$ score of the temporary candidate network no longer decreases, indicating convergence. The ${\tilde{G}}_k$ after the above procedure is regarded as the final potential optimal network in this repetition.

Step 4. Determine a global optimal network.

After obtaining potential optimal networks from different repetitions $\left\{{\tilde{G}}_1,\dots ,{\tilde{G}}_I\right\}$, we take the network with the minimum $H$ score among them as the global optimal network $\tilde{G}$, that is:

(14) $$\begin{array}{c}\tilde{G}=\underset{{\tilde{G}}_k\in \left\{{\tilde{G}}_1,\dots ,{\tilde{G}}_I\right\}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\{H({\tilde{G}}_k)\},\end{array}$$thereby qualitatively and quantitatively describing the system’s network state at the current time point.

Step 5. Identify critical states based on a dual evaluation mechanism of network scoring and structural stability.

At each time point $t$, we calculate the $H$ score of the global optimal network $\tilde{G}\left(t\right)$ (i.e., $H\left(\tilde{G}\left(t\right)\right)$) following the above procedures. This score, named BCTI score and denoted by $S\left(t\right)$, provides a quantitative measure of the network’s temporal dynamics and serves as a key indicator for identifying potential critical transitions. In terms of the DNB theory, the system approaching a critical point typically exhibits distinct dynamical changes within a core subnetwork: (i) sharply increased fluctuations (such as standard deviation) in the expression of key genes; (ii) markedly enhanced coordinated variation (such as Pearson correlation coefficient) among key genes. In our method, these two features are reflected in the dimensions of network scoring and structural stability, respectively, and jointly drive the identification of critical states.

Specifically, in the dimension of network scoring, it is demonstrated that the defined scoring function $H(\cdot )$ is theoretically derived to be positively correlated with the variance of gene expression (Note S3). Due to the markedly increased gene expression fluctuations near the critical state, the BCTI score rises correspondingly and can therefore serve as a quantitative indicator of declining system stability. To statistically identify the effective signal, we apply a one-sample t-test to the BCTI scores. The time point $T=t$ is defined as a critical point if there is a significant difference between the current BCTI score $S\left(t\right)$ and the mean value of a vector $\left(S\left(1\right),S\left(2\right),\dots ,S\left(t-1\right)\right)$ ( $P<0.05$). More details are provided in Note S4. In the dimension of structural stability, BCTI employs structural equation models to reconstruct gene regulatory networks. When the system is at the critical point, markedly enhanced coordinated variation among key genes leads to severe multicollinearity, which substantially increases the variance of regression coefficient estimates during structural equation model fitting (Grewal, Cote & Baumgartner, 2004). It disrupts the stable inference of causal relationships and leads to instability in network reconstruction, thereby reducing the accuracy of identifying regulatory relationships. Therefore, the combination of elevated BCTI scores and decreased network stability provides a dual criterion for identifying an impending critical transition and helps distinguish true system transitions from isolated anomalies.

Efficiency of BTCI on inferring causal relationships in GRNs

To evaluate the performance of BCTI in inferring causal relationships within GRNs, we compared it against several methods, including GENIE3-RF-sqrt, GENIE3-ET-sqrt, GENIMS, GNIPLR, NARROMI, NIMEFI, and PLSET, using benchmark datasets such as the DREAM network, the IRMA network, and the SOS DNA repair experimental dataset. For the DREAM and IRMA networks, BCTI successfully inferred most of the edges in the gold standard networks, achieving 85.56% and 85.00% accuracy, respectively, outperforming other methods (Tables S3 and S4). On the SOS DNA repair dataset, BCTI correctly inferred 14 out of 24 regulatory edges and achieved the highest true positive (TP) count, indicating strong performance in recovering the underlying regulatory structure (Table S5). These results demonstrate the effectiveness and robustness of BCTI in inferring gene regulatory networks across both simulated and experimental datasets, highlighting its ability to achieve accurate and reliable causal relationship inference under diverse data conditions. Additional computational details are provided in Tables S3–S5.

Accurate and robust critical-state detection performance of BCTI compared with existing methods

To examine the ability of BCTI to identify critical transitions, we employed a simulation framework based on a previously published gene regulatory network model that has been widely used to study bifurcation-driven dynamics in biological systems (Foo, Kim & Bates, 2018; Chen et al., 2015). In this framework, regulatory interactions are described using Michaelis–Menten–type kinetics and the system evolution is governed by stochastic differential equations. Following this established simulation setting, system dynamics were generated by gradually varying a control parameter ( $q$) within the range of −0.3 to 0.3, which drives the system toward a bifurcation point. The resulting simulated time-series data were then used to evaluate the effectiveness of BCTI in detecting the onset of critical states. Detailed model formulations and parameter settings are provided in Note S6.

As depicted in Fig. 2B, the BCTI score exhibited an abrupt increase as the system approached to the bifurcation point ( $q=0$), indicating an impending critical point. To visually emphasize the dynamic differences between the before-transition and critical states, we presented the evolution of the network reconstructed by BCTI (Fig. 2C). Notably, when the system was far from the critical point (such as $q=-0.2$), the reconstructed network exhibited notable overlap ( $\mathrm{T}\mathrm{P}=12$) with the gold standard network, highlighting the system’s stable state and demonstrating BCTI’s efficacy at such a stage (Fig. 2D). However, as the system was near the critical point ( $q=-0.001$), there was a notable decrease in the overlap ( $\mathrm{T}\mathrm{P}=6$) between the reconstructed network and the gold standard network, indicating an impending transition of the system’s state (Fig. 2D). To further underscore the robustness of BCTI, we conducted comparative experiments under varying noise levels against several benchmarking methods (Chen et al., 2012; Peng et al., 2022), including the node-based, direct interaction network-based divergence (DIND), and DNB methods (Fig. S5). Under noise-free conditions (noise strength $\sigma =0$), all methods were able to indicate the tipping point. However, as the noise strength increased (e.g., $\sigma =0.5$ or $\sigma =1.5$), the node-based, DIND, and DNB methods exhibited substantial instability or even failed to identify the critical transition effectively. In contrast, BCTI consistently exhibited its effectiveness and robustness in identifying tipping points, maintaining high sensitivity and producing clear early-warning signals across all tested noise levels. More details for the theoretical background of BCTI and computational results are provided in Note S3 and Table S6. The above numerical experiment validates the reliability and accuracy of BCTI in identifying critical states at the network level.

Efficiency of BCTI on uncovering underlying gene regulatory mechanisms of cancers

Furthermore, BCTI was applied to three TCGA datasets: COAD, LUAD, and THCA. For each cancer type, we selected biologically relevant pathogenic pathways known to be strongly associated with the disease and extracted key genes from these pathways to reconstruct the gene regulatory networks. In other words, these genes were predefined prior to network inference based on curated biological knowledge. The same set of genes was used across all tumor stages within each cancer type to ensure the comparability of the reconstructed networks. For COAD, the verification network was derived from the pathway hsa04151 (PI3K-Akt signaling pathway) in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.genome.jp/kegg/pathway.html). For LUAD, we took the pathway hsa04064 (NF-kappa B signaling pathway) as the verification network due to its tumor promoting role in lung carcinogenesis (Chen et al., 2011). For THCA, the pathway hsa04010 (MAPK signaling pathway) was selected as the verification network, as it has been widely studied and plays a crucial role in thyroid carcinogenesis (Hu et al., 2021). Implementing the procedure outlined in the Methods, we rewired the genetic network within the functional pathways at each tumor stage and calculated BCTI scores to identify the critical states during cancer progression (Fig. 3). On one hand, as shown in Figs. S6–S8, BCTI characterized the dynamic evolution of the rewired regulatory subnetworks in the PI3K-Akt signaling pathway across all five stages, the NF-kappa B signaling pathway across six stages, and the MAPK signaling pathway across four stages, respectively. Overall, BCTI demonstrates strong network reconstruction capabilities by accurately identifying true regulatory edges within pathways during early stages of disease progression. However, its performance declines markedly at certain critical stages. For instance, only two true regulatory edges were detected in stage IIB of COAD ( $TPR=0.1333$), none were identified in stage IIIB of LUAD ( $TPR=0$), and only three in stage II of THCA ( $TPR=0.1250$). Compared to other stages, the notable drop in reconstruction accuracy in these stages suggests the system approaches a critical state. According to the network stability criterion proposed in the Methods section, this phenomenon may arise from enhanced coordinated variation among key genes near the critical state, which increases multicollinearity among variables and consequently disrupts the stable estimation of regulatory relationships in the structural equation models. More computational details are provided in Table 1.

On the other hand, the significant increase in the BCTI score identified the critical state as stage IIB for COAD ( $P=0.0346$), stage IIIB for LUAD ( $P=0.0006$), and stage II for THCA ( $P=0.0302$) (Fig. 3A). These results were consistent with those reported in previous studies (Liu et al., 2022; Zhong et al., 2023). However, the gray curves derived from the mean expression of differentially expressed genes (DEGs) and the genes involved in the networks demonstrated that they cannot signal such critical transitions (more computational details are provided in Note S7 and Fig. S9). Moreover, we also conducted comparative experiments to evaluate the performance of BCTI on the real datasets against other representative approaches, including the edge-based relative entropy (ERE), DNB, and DIND methods (Chen et al., 2012; Peng et al., 2022; Hong et al., 2024) (Figs. 3B–3D and Table S7). The results suggested that compared to other methods, the critical points inferred by BCTI exhibited higher statistical significance (lower $P$-values) and a closer correspondence to clinical observations. Specifically, for COAD, cancer metastasis has not occurred at the identified critical state (stage IIB); however, by stage III, tumor cells have already reached the nearby lymph nodes (Amin et al., 2017). For LUAD, cancer has not metastasized at the identified critical state (stage IIIB), whereas it has already spread to distant tissues or organs via the bloodstream at stage IV (Chiang & Massagué, 2008). For THCA, cancer metastasis is absent at the identified critical state (stage II); however, by stage III, cancer has spread to regional lymph nodes (Shaha, 2007). Whereas the ERE, DNB, and DIND methods failed to capture such statistically and biologically meaningful critical points.

To validate the identified critical state, we performed Kaplan–Meier (log-rank) survival analysis for samples separately from before and after the identified critical stage. As depicted in Fig. 3E, the survival expectancy is much higher before the identified critical state than after it, with significant $P$ values observed for COAD, LUAD, and THCA ( $P<0.0001$, $P=0.0019$, and $P<0.0001$, respectively). Moreover, the prognosis analysis also supports the computational results based on BCTI. For example, compared with the prognosis analysis based on other stage divisions, the difference in survival expectancy before and after the identified critical state of THCA by BCTI (stage II) is much more significant ( $P<0.0001$) (see Fig. S10 for details of the prognosis analysis). These results highlight that patients diagnosed before the identified critical state by BCTI have significantly better prognoses than those diagnosed after the critical stage. More detailed information regarding the survival analysis of tumors can be found in Note S8 and Figs. S11, S12. Taken together, the BCTI score can identify early warning signals for critical states associated with survival expectancy.

In addition, for THCA, we performed Gene Ontology (GO) functional enrichment analysis on differentially expressed genes at each stage. It was found that some regulatory and target genes (such as IL1A and TNF) (Figs. 3F), although not included in the gold standard, were simultaneously significantly enriched in the same tumor-associated GO terms (such as negative regulation of apoptotic signaling pathway (GO:2001234)) (Figs. 3G and Fig. S13), suggesting the potential biological significance of the regulatory relationships inferred by BCTI (see Tables S8 and S9 for more details). Indeed, in the early stage of THCA, IL1A may induce the expression of TNF, which could trigger further cascades of inflammatory factors (Guarino et al., 2010). This process might promote tumor cell proliferation and enhance their anti-apoptotic capabilities. Moreover, IRAK1 and MYD88 were simultaneously significantly enriched in the regulation of I-kappaB kinase/NF-kappaB signaling (GO:0043122) (Fig. 3F and Table S8). Biologically, chronic inflammation-driven IRAK1 activity may upregulate MYD88 expression through pathways such as NF-kappaB, thereby amplifying inflammatory cascades and promoting tumor progression (Jain, Kaczanowska & Davila, 2014). All the results revealed insights into the microscopic change of the underlying signaling transition mechanism, demonstrating BCTI’s capability to offer new perspectives for understanding cancer progression mechanisms and developing precision medicine strategies.

Leveraging BCTI on single-cell datasets to uncover gene regulatory mechanisms in lung development

For the case study of BCTI on human embryonic lung development data, we utilized single-cell data from the pseudoglandular period at 8, 10, 11.5, 12, 13, and 14 post-conception week (PCW), including a total of 42,982 cells and seven main cell types (including mesenchymal, megakaryocyte, epithelial, endothelial, immune, neuronal, and erythroblast/red blood cells) (Fig. 4A) (Sountoulidis et al., 2023).

First, highly variable genes (HVGs) were identified for each stage separately, and their union was used to construct the gene regulatory network. The BCTI method was then applied to quantitatively analyze the dynamic changes in the network, thereby identifying the critical state during lung development (Fig. 4B). As shown in Fig. 4C, the significant increase in the BCTI score identified the critical state as 10 PCW ( $P=0.0256$). However, the mean expression of the HVGs failed to provide the effective early-warning signal of critical transition during lung development (the gray curve).

Subsequently, key hub genes, ranked by degree, were selected based on degree centrality analysis in the networks before and after the critical state, and biological functional analysis was performed on these genes to explore their downstream functional impacts caused by major regulatory transitions (Fig. 4B). The results showed that, on one hand, in terms of network structure, compared to the before-transition-state GRN, the after-transition-state GRN exhibited tighter connectivity among key hub genes and a higher degree centrality of these genes (Fig. 4D and Fig. S14). It suggested a more consolidated and cooperative regulatory architecture, where key hub genes played an increasingly central role in network stability and functional coordination. On the other hand, in terms of biological function, key hub genes before the critical state were predominantly enriched in functional modules associated with epithelium-related lung development, such as lung epithelium development (GO:0060428) and lung epithelial cell differentiation (GO:0060487), facilitating lung morphogenesis and functional differentiation (Kimura & Deutsch, 2007) (Fig. 4E). However, after the critical state, key hub genes were primarily enriched in functional modules related to endothelium-related lung development, such as morphogenesis of a branching structure (GO:0001763) and endothelial cell development (GO:0001885), highlighting the growing importance of endothelial signaling in supporting lung alveolarization, angiogenesis, and the establishment of the alveolar-capillary network, which are essential for functional gas exchange (Lignelli et al., 2019) (Fig. 4E). Additionally, the expression patterns of certain key hub genes exhibited cell-type heterogeneity, such as GPC3 and SOX2 at 8 PCW and CCBE1 and CDH5 at 11.5 PCW, further supporting the findings of the GO functional analysis (Fig. 4F). Indeed, the above results aligned with the previous study showing that foetal breathing movements, which are considered to play important roles in lung development (Nikolić, Sun & Rawlins, 2018), begin at 10-11 PCW. More analytical details of BCTI are provided in Table S10. These findings underscore the ability of BCTI to capture dynamic changes in gene regulatory networks that correspond to critical developmental processes.