Computational strategies in systems-level stress response data analysis

1.1 Cellular homeostasis and stress response

A principal challenge in stress research is establishing a concise, universal, and unambiguous definition of “stress”. Conceptually, however, stress in biological systems can be reduced to a two-state transition model (Kültz 2020). In this model, stress is a force acting upon an organism that shifts it from a “normal” homeostatic state to a stress-induced state of dysregulation (Figure 1). Under standard or “normal” conditions – representing an evolutionary average of environmental and genetic forces – cells and organisms exist in a dynamic equilibrium. This equilibrium can be maintained through the energy-dependent cellular homeostasis response (CHR), which counteracts deviations in biological set-points arising from continuously changing environmental conditions (Chovatiya and Medzhitov 2014; Galluzzi et al. 2018; Somero 2020).

Figure 1:

Simplified transition state model of cellular responses to stress. Under normal conditions, the cellular homeostasis response maintains a biological system in stable equilibrium or homeostatic state. When stress forces surpass the system’s buffering capacity, the system transitions into a stress state. During this phase, the cellular stress response activates mechanisms to counteract and mitigate strain. If successful, the system can revert to homeostasis, where the cellular homeostasis response again preserves stability; otherwise, permanent damage may lead to irreversible cell death.

Nevertheless, every biological system possesses a finite buffering capacity, or elasticity that keeps physiological parameters, such as rates, within homeostatic limits. When stress exceeds this buffering capacity, a transition to a cellular stress response (CSR) occurs. Like the CHR, the CSR is also energy dependent, but it is specifically triggered when (macro)molecular integrity is compromised, forcing the system out of its homeostatic range. In this state, macromolecular damage can intensify overall dysregulation while further activating CSR pathways (Galluzzi et al. 2018; Simmons et al. 2009). For instance, molecular sliding clamps monitor DNA damage and initiate repair upon surpassing a critical threshold, whereas heat shock proteins activate transcription factors that address protein unfolding (Freeman and Monteiro 2010; Somero 2020). If the combined efforts of the CHR and CSR fail to restore the system to homeostasis, the resulting progressive imbalance may culminate in system failure.

1.2 Stress effect and response specificity in biological systems

Stress, defined as the sum of forces acting on a biological system, triggers a range of responses that enable the organism to counteract and mitigate strain. Due to the large number of variables that can be affected by stress, spanning genes, proteins, metabolites, and their interactions, stress responses often exhibit remarkable complexity and context-dependence while changing dynamically over time (Boos et al. 2020; Pommerrenig et al. 2018; Schroda et al. 2015; Sies et al. 2017). Moreover, the magnitude of the stressor and the duration of exposure converge to determine whether a system remains within its elastic limit or transitions irreversibly into a stress state. Therefore, stress responses also differ substantially depending on the type of stress (Garcia-Molina et al. 2020; Kreps et al. 2002). Eukaryotic cells possess an impressive capacity to cope with diverse environmental challenges (Boos et al. 2019; McShane et al. 2016; Ronen and Botstein 2006). Rapid and often reversible adjustments to chemical or osmotic changes, temperature fluctuations, variations in light intensity, shifts in metabolic demands, and nutrient availability are mediated by intracellular response programs (Fulda et al. 2010; Karpinski et al. 1999; Kourtis and Tavernarakis 2011; Ronen and Botstein 2006). While many of these programs share fundamental regulatory principles, for example, upregulation of chaperones, ubiquitin-proteasome machinery, redox enzymes, and global adjustments in gene expression, their understanding of how different stressors activate distinct subsets of these routines or recruit entirely additional, specialized pathways is challenging (Kourtis and Tavernarakis 2011; Pommerrenig et al. 2018; Schroda et al. 2015; Sies et al. 2017). Consequently, stress research demands a context-specific perspective, requiring scientists to systematically identify the molecules, pathways, and interactions most critical in each situation.

Elucidating stress responses at a system-wide level relies on two key steps: (1) defining the modules and components involved – along with their identity, quantity, and timing of action – and (2) analyzing how these components interact within a molecular network to characterize responses and uncover key regulatory links. By integrating empirical observations with theoretical models, researchers progressively reveal how biological systems maintain homeostasis under various stress intensities (Figure 2). This knowledge illuminates fundamental physiological principles and supports the development of interventions targeting the most vulnerable aspects of stress responses. Consequently, computational systems biology approaches and especially network modeling are indispensable for mapping the underlying “wiring diagrams” of cellular regulation and pinpointing pivotal nodes and interactions.

Figure 2:

Analysis of stress effects and response specificity. Stress responses vary significantly depending on the stressor type, activating distinct intracellular subroutines that ensure specificity. Computational models help define the modules and components involved, including their identity, quantity, and timing, which clarifies the dynamic nature of these responses and reveals key regulatory mechanisms.

2.1 Statistical analysis and testing

When a single stress time point is of primary interest, a classical two-sample test (such as a t-test or its non-parametric equivalent) can determine whether the mean value at that time point diverges significantly from control levels (Gosset 1908; Helmert 1876; Pearson 1895). In many contemporary study designs, molecule abundances are compared within a single sample over time or between treated and control samples. Thus, one can interpret a molecule’s time series as a continuous shift from control to treatment and potentially back again (Garcia-Molina et al. 2020; Zhang et al. 2022). For each protein, transcript, or other molecule of interest, a statistical procedure tests whether it changes significantly at any measured point.

A widely used method for detecting overall changes across multiple time points is one-way Analysis of Variance (ANOVA) (Fisher 1992; Stahle and Wold 1989). This test evaluates data from all replicates at each time point to identify significant differences in group means. A more refined option, repeated measures ANOVA, accounts for correlations among measurements taken from the same sample, but relies on assumptions (e.g., sphericity) that may not hold in extended studies with missing or imputed data (Keselman et al. 2001; Park et al. 2009). Consequently, some researchers opt for simpler frameworks – such as the standard one-way ANOVA – though these often have reduced sensitivity to time-dependent effects (Knöringer et al. 2023; Mergner et al. 2020; Sacco et al. 2016). Post hoc tests after ANOVA can indicate which specific groups or time points differ significantly from each other (Armstrong 2014; Agbangba et al. 2024; Brown 2005). In stress research, they are often used to infer the sequence of events, as demonstrated through the identification of targets for enhancing plant resistance to abiotic stress in Oryza sativa, revealing differences in susceptibility to mitoprotein-induced stress in yeast, allowing the characterization of proteome-wide changes at different diabetic stages in mice or by the investigation of cell type-specific alterations in ABA signalling pathways in response to alterations in salinity (Geng et al. 2013; Knöringer et al. 2023; Moin et al. 2016; Sacco et al. 2016). While powerful, standard statistical approaches cannot easily distinguish between the amplitude and timing of changes. Moreover, many such tests assume homoscedasticity (equal variances across time points), making them unsuitable for highly variable datasets. While non-parametric options exist, having only three or four replicates can weaken statistical power and inflate false positives. Consequently, researchers must weigh theoretical assumptions against practical constraints. Depending on data design and goals, methods like Dunnett’s test, repeated measures ANOVA variants, or mixed-effect models offer more robust outcomes (Dunnett 1955; Laird and Ware 1982; Park et al. 2009; Wood 2013). However, these approaches do not fully account for the strong dependence and ordering inherent in time series data, highlighting the need for more specialized models that explicitly address temporal dependencies (Jung and Tremayne 2003; Kin Kwan Leung et al. 2023).

2.2 Model-based time series representation

As high-throughput technologies become more affordable, researchers can increasingly incorporate the temporal dimension into stress response studies. Unlike datasets without inherent ordering, time series data exhibit strong dependencies among consecutive measurements, with each data point influenced by its predecessors and providing context for subsequent observations (Jung and Tremayne 2003; Kin Kwan Leung et al. 2023). In principle, if environmental conditions remain unaltered, one might expect values such as molecule abundances to fall within a range bounded by neighboring time points. However, this assumption proves overly simplistic in living systems, especially under stress, because internal regulations, ongoing processes, and expression noise can induce fluctuations that are neither linear nor easily predictable (Cohen et al. 2008; Elowitz et al. 2002).

Biological noise complicates modeling efforts because it cannot be fully controlled or anticipated (Chowdhury et al. 2021). Regulatory cascades often depend on threshold effects, making sudden “bursts” in gene expression or protein levels possible (Beckman et al. 2021; Gardner et al. 2000; Luo et al. 2023). Predicting such events requires complete knowledge of every factor influencing a molecule’s state. While great strides have been made in optimizing analysis and prediction algorithms, often by focusing on smaller, well-defined processes, fully modeling the multitude of molecules involved in an entire stress response remains challenging. For instance, metabolic flux models can accurately predict changes in metabolite concentrations by tracking enzyme activities and co-factor abundances in isolated pathways (Kim et al. 2008; Weaver et al. 2014). Yet, applying such detailed mechanistic models to large-scale stress responses is limited by missing information on protein-coding genes, unknown interactions among nucleic acids, proteins, and metabolites, and variability introduced by genetic heterogeneity and imperfect measurement (Fröhlich et al. 2018; Ghatak et al. 2019; Lee et al. 2012).

Consequently, more general model-based time series representations are often employed in stress response research. Among these, clustering techniques remain popular because they impose relatively few assumptions about the underlying biological processes. Methods such as k-means clustering of time-point averages, as well as hierarchical or density-based clustering, have successfully revealed dynamic differences within and across responses (Datta and Datta 2006; Maigné et al. 2023; Tavazoie et al. 1999; Wang et al. 2008). By grouping similar temporal trajectories, these approaches provide insight into how cells coordinate molecular changes under varying stress conditions, even when complete mechanistic details are unavailable. This is demonstrated in studies of Medicago truncatula under drought stress, where k-means clustering of metabolite concentrations identified tissue-specific responses, notably in sucrose and citrate levels in the leaves (Dickinson et al. 2018). Similarly, k-means clustering to analyze both gene expression and protein abundance patterns during Solanum lycopersicum pollen development under heat stress revealed stage-specific heat stress responses, with tetrad-stage pollen showing the strongest response and mature pollen displaying minimal changes (Keller and Simm 2018).

However, traditional clustering methods usually treat each time point as independent, which can obscure crucial temporal dependencies. For instance, time series vectors might appear equally distant even if their order is altered – a scenario at odds with biological intuition, where each measurement often depends on previous states (Abanda et al. 2019; Nies et al. 2019). To address these limitations, distance measures that consider temporal context, such as Dynamic Time Warping (DTW), have been introduced (Berndt and Clifford 1994; Müller 2007). DTW4Omics, a tool based on DTW, identifies associations between biological entities and endpoints by dynamically aligning their time series, even when their temporal patterns are not synchronized. When applied to existing data from a study on oxidative stress effects in human colon carcinoma cells, DTW4Omics revealed novel gene associations that conventional correlation analyses missed (Cavill et al. 2013). Newer approaches aim to represent individual molecular trajectories, using interpolation or regression to incorporate information from neighboring time points. Constrained smoothing splines, for example, enforce biological plausibility (e.g., by limiting excessive oscillations) and can correct for time-dependent variance (Venn et al. 2024). Collectively, these methods enhance the interpretation of temporal data and refine our understanding of complex regulatory mechanisms.

2.3 Identification of stress response modules via enrichment analysis

Time-series modeling and statistical analyses offer a concise view of stress response data. In parallel, the availability of functional or physicochemical information about genes and molecules can help define specific response groups more accurately. Leveraging these complementary data sources allows researchers to pinpoint critical “response modules” and determine how particular pathways become activated or suppressed during stress.

A central approach to interpreting lists of differentially regulated molecules is Gene Set Enrichment Analysis (GSEA) (Subramanian et al. 2005). GSEA identifies sets of molecules overrepresented in functionally coherent pathways, compared to a broader background of measured entities. Most often, these sets derive from established annotation databases such as MapMan, GO, KEGG, Reactome, Wikipathways, or BioCyc (Ashburner et al. 2000; Jassal et al. 2020; Kanehisa and Goto 2000; Karp et al. 2005; Kelder et al. 2012; Thimm et al. 2004). A common method for performing GSEA is the one-sided hypergeometric or Fisher’s exact test, which detects significantly enriched sets in an experiment (Al-Shahrour et al. 2004; Da Huang et al. 2009; Maere et al. 2005; Zeeberg et al. 2003; Zhou and Su 2007; Zhang et al. 2004; Zhong et al. 2004). Here, each molecule is assigned a p-value or label indicating a significant change over time or relative to a control, followed by a hypergeometric test to identify enriched sets (Rivals et al. 2007). Because each term necessitates a separate test, multiple-testing corrections are crucial. GSEA is frequently applied in the context of stress research, as numerous examples illustrate. In yeast, GSEA allowed two distinct transcription factor groups which induce sharply different gene sets, reflecting two distinct cellular strategies for coping with escalating levels of H₂O₂ stress. In human neutrophils, GSEA revealed that stimulus-specific transcriptional shifts (for example, G-CSF versus IFN-driven) intertwine with neutrophil maturation stage, leading to clearly different gene-expression signatures under stress myelopoiesis. In O. sativa it allowed to pinpoint that wet conditions favor the selection of growth and defense genes, while drought conditions lead to an enrichment in genes associated to the water stress response, growth adaptation, and flowering regulation (Groen et al. 2020; Jose et al. 2024; Montaldo et al. 2022). However, relying solely on p-value thresholds leads to a binary categorization that can obscure nuanced information, a problem intensified in time-series data (Pan et al. 2005). To address this, alternative GSEA methods are often divided into two main categories. Functional Class Scoring (FCS) aggregates ranks or p-values of individual molecules within a set, providing an overall significance score for that set. This strategy helps avoid arbitrary cutoffs by assessing whether the molecules in a pathway collectively shift in response to stress (Pavlidis et al. 2002; Tarca et al. 2013). In contrast, Single-Sample (SS) methods compute an enrichment score for each pathway in each sample, highlighting which pathways are most relevant to specific samples or time points (Barbie et al. 2009; Dinu et al. 2007; Frost et al. 2015; Shen and West 2008; Tarca et al. 2013). By assessing gene set activity on a per-sample basis, SS methods can reveal subtle, context-dependent variations in pathway regulation that might remain hidden when data are pooled. In stress research, the integration of SS GSEA with single-cell sequencing technology was employed to investigate the differential response to oxidative stress between gastric cancer and normal gastric tissue. By leveraging enrichment scores calculated by SS GSEA from GO biological process terms in the Molecular Signature Database, comparative analyses revealed that gastric cancer tissue, exhibits diminished oxidative stress response capabilities. Succeeding examination of the gastric microenvironment at the single-cell level demonstrated cell type-specific responses, particularly in T-cells (Yu et al. 2022).

Further refinement of GSEA involves integrating multiple annotation databases or addressing the issue of overlapping gene sets, acknowledging that the same molecule can participate in different biological pathways and processes (Simillion et al. 2017). These overlaps can introduce analytical complexities (e.g. double-counting or misattribution of gene function) which, if not properly addressed, may obscure the true extent of pathway activation. Methods that merge or reconcile annotations from diverse databases thus provide a more robust view of how stress responses are coordinated at the systems level.

Enrichment analyses are primarily statistical, but accounting for biological realities can enhance understanding of stress responses. This is especially relevant in scenarios where a small number of pivotal enzymes exert disproportionate control over entire pathways. Biophysically motivated models, which estimate the “investment” in each pathway, can illuminate the most crucial factors driving stress responses (Schneider et al. 2020). By focusing on these pivotal elements, researchers gain deeper insights into how specific pathways contribute to stress tolerance and resilience.

3.1 De novo inference of stress response networks

De novo inference begins without a predefined network structure, often relying on large-scale datasets, published literature knowledge or node features derived from computational predictions, to assemble networks (Luck et al. 2020; Szklarczyk et al. 2023). Researchers employ various strategies to construct stress response networks de novo. These approaches range from pairwise methods, which focus on relationships between node pairs, to context enriched methods that leverage overall network topology, e.g. incorporating neighborhood structures, paths, and other contextual information to infer more accurate and robust networks (Figure 4).

Figure 4:

Workflow for de novo inference of stress-response networks. (i) Multi-omics measurement under distinct stress conditions produces quantitative read-outs such as RNA-seq, proteomics, and metabolomics. (ii) These raw signals are transformed into a condition-specific feature matrix (heat map), optionally augmented with contextual evidence from curated pathways or literature mining. (iii) Pairwise and topology-aware algorithms then generate an ensemble of candidate networks. (iv) Finally, statistical and biological criteria are applied to this ensemble to select a single, parsimonious stress-response network model.

Pairwise approaches for de novo network construction involve the application of measures like correlation or mutual information measures to sets of individual nodes. Oftentimes each node (e.g., corresponding to a gene or protein identifier) is associated with a set of numerical values (a.k.a. features) derived from omics datasets e.g. a genes abundance over a time course or computational predictions, sequence or structure similarity. In their simplest form, these networks are constructed by relating individual nodes through pairwise measures and linking only those where the metric exceeds a chosen threshold, as in WGCNA (Langfelder and Horvath 2008). However, more sophisticated methods for threshold determination exist, incorporating global features or neighborhood information. Notable examples include ARACNe and CLR (Faith et al. 2007; Margolin et al. 2006). The resulting networks are typically undirected graphs. Whether these networks are signed or unsigned depends on the specific metric and thresholding applied. For instance, Pearson correlation is frequently used to sign the edges, whereas in some cases, the absolute value is computed, resulting in undirected edges. In the context of stress research, the application of WGCNA to expression profiles derived from Chlamydomonas reinhardtii cells responding to 35 °C and 40 °C exposures over time allowed to reconstruct the differential regulation of key biological pathways with respect to different temperatures. The resulting adjacency matrix was used to identify modules, which led to the identification of twelve key transcriptomic modules, each showing characteristic temporal expression patterns during heat treatment and recovery (Zhang et al. 2022).

Additionally, node features derived from conditional omics data can be enriched to introduce causality by integrating evidence such as chromatin immunoprecipitation sequencing (ChIP-Seq), protein interaction assays, literature-derived annotations, or computationally predicted binding sites (Chow et al. 2019). For example, distinct protein classes, such as transcription factors versus target genes, can be represented as separate node types in heterogeneous networks, or incorporated into the node feature space to capture causal dependencies (Wang et al. 2022b). In doing so, edges can be modelled as directed relationships, as demonstrated in the simultaneous use of ChIP-Seq and differential expression data to identify context-specific regulatory relationships, such as those modulated by environmental stresses. For example, in the context of the abscisic acid (ABA) signaling pathway in Arabidopsis thaliana, genome-wide ChIP-Seq profiling of ABA-responsive transcription factors (TFs) was integrated with DREM (Dynamic Regulatory Events Miner) analysis of time-series RNA-seq data to construct a comprehensive, time-resolved ABA-specific gene–TF network. TFs were subsequently arranged into hierarchical tiers based on their connectivity (in-degree and out-degree), illuminating which were higher-level regulators. This approach uncovered novel ABA-related genes and underscored how coordinated, dynamic binding events and hierarchical regulation jointly shape the plant’s stress-responsive transcriptional program (Song et al. 2016a).

Since de novo networks are often noisy and may miss key interactions or hidden relationships, refinement steps such as link prediction and node classification are essential to improve network accuracy. These steps also help reveal emergent properties, such as modules or communities, which can provide deeper insights into functional organization and regulatory mechanisms within stress response networks.

3.2 Refinement of stress response networks

Network refinement encompasses a broad set of tasks such as classifying or grouping nodes and modifying edges based on new experimental or computational insights, by drawing on contextual graph information and additional biological data (Liu et al. 2020). These efforts include annotating previously uncharacterized nodes, assigning nodes to modules, reweighting or removing links, and predicting missing links (Figure 5). Moreover, refinement often accounts for dynamic or context-specific influences, enabling networks to capture the fluctuations and adaptive responses that biological systems exhibit under stress. Classical network refinement often relies on topological algorithms to detect communities or modules within the network. Community detection algorithms, such as the Girvan-Newman and the Louvain algorithm, identify clusters of nodes that are more densely connected internally than with the rest of the network (Blondel et al. 2008; Girvan and Newman 2002). These methods leverage the inherent structure of the network to uncover functional modules, which can correspond to biological pathways, protein complexes, or regulatory circuits. In the context of stress response research, the application of the Louvaine method allowed the identification of distinct regulatory states within one cell population, by the refinement of an de novo generated co expression network based on a shared nearest-neighbor graph (Jackson et al. 2020).

Figure 5:

Refinement operations applied to an initially inferred stress-response network. Starting from the raw graph, (i) latent-space link prediction evaluates node-embedding vectors (v₁, v₂); when their similarity exceeds a learned threshold, a new edge is introduced, and (ii) neighborhood-based grouping then assigns the hub’s adjacent nodes to putative functional modules, indicated by arrows drawn parallel to the edges. The right-hand panels display the results: the upper panel highlights newly added, high-confidence links that enrich the topology, whereas the lower panel shows how module assignment delineates previously uncharacterized communities.

Today, machine learning (e.g. graph learning) approaches have augmented classical methods by learning graph embeddings based on network topology as well as integrating node features into the refinement process. Graph embedding methods, such as node2vec and DeepWalk, use topological information to transform the network into a latent space where the proximity of node embeddings reflects their likelihood of interaction (Grover and Leskovec 2016; Perozzi et al. 2014). Graph Neural Networks (GNNs) can incorporate diverse biological attributes encoded as node features, including gene expression patterns, functional annotations, and domain-specific knowledge and propagate this feature information across the network, allowing the model to learn representations that capture both the local topology and the intrinsic properties of the nodes (Depuydt and Vandepoele 2021; Hetzel et al. 2021; Wang et al. 2022a). Thus, when experimental data can be associated with specific conditions, refinement methods can incorporate this conditional information to context-specific node classifications as demonstrated by GeneWalk, which associates genes to functions relevant in the biological context under study (Ietswaart et al. 2021). GeneWalk demonstrates this capability in the analysis of transcriptional responses to the bromodomain inhibitor JQ1 in T-cell acute lymphoblastic leukemia cells. While traditional GO enrichment analysis only identified five generic functions like “ncRNA metabolic process” and “chromatin organization” with low fold enrichment, GeneWalk’s network-based approach revealed condition-specific functions for individual genes. For example, when examining BRCA1, which was downregulated after JQ1 treatment, GeneWalk ranked DNA damage and repair-related processes as most relevant, followed by histone modifications, while other annotated functions like transcription and metabolism were deemed less significant in this context (Ietswaart et al. 2021).

Equally crucial to node refinement is the refinement of network edges. This can involve reweighting existing links based on new evidence or removing spurious or low-confidence interactions. Regarding co-expression networks methods which build on Random Matrix Theory (RMT) have been used to distinguish signal from noise without including information of a networks structure a priori in conditional networks (Ficklin et al. 2017; Gibson et al. 2013; Luo et al. 2006).

Likewise, embeddings derived by machine learning approaches can be used to predict missing links by identifying node pairs with high similarity in the latent space but lacking an existing edge (Chen and Liu 2022; Guo and Xiao 2024; Mao et al. 2023). Graph neuronal network-based approaches like GENELink or GNNLink utilize Graph Attention Networks or Graph Convolutional Neuronal Networks to integrate prior knowledge of gene interactions with gene expression data derived from single-cell RNA sequencing (Chen and Liu 2022; Mao et al. 2023). Moreover, Directed Graph Convolutional neural network-based method for GRN inference (DGCGRN) further integrates genes sequence features to iteratively update links in an given gene regulatory network, and evaluates its capabilities on Escherichia coli data sets subject to cold-, heat- and oxidative-stress (Wei et al. 2024). A case study applying DGCGRN to reconstruct the bladder urothelial carcinoma network displayed its effectiveness, DGCGRN outperformed algorithms like GENIE3 and PLSNET, while also identifying novel regulatory interactions that were validated through literature research (Wei et al. 2024).

3.3 Comparative analysis of stress response networks

Comparative network analysis provides a powerful lens for characterizing condition-specific stress responses by examining and integrating multiple biological networks to identify shared or unique components. A range of computational approaches can be used to detect rewiring events and conserved substructures across networks derived under different conditions (Emmert-Streib et al. 2016; Tantardini et al. 2019).

When node sets are shared across the networks, such as condition-specific networks generated for the same species, global (e.g. degree distribution, modularity) and local (e.g. node centrality) network statistics can be used for direct comparisons (Figure 6). A direct comparison can be achieved by formulating proximity measures which make differences and similarities between networks quantifiable (Manipur et al. 2020). Besides the straightforward application of measures like the Pearson correlation coefficient Euclidean, Manhattan, or Jaccard distance to adjacency matrices, proximity measures like DeltaCon and the CutDistance have been developed (Emmert-Streib et al. 2016; Koutra et al. 2013; Liu et al. 2018; Tantardini et al. 2019; Tesson et al. 2010). A combination of both approaches can be found studying A. thaliana, where condition-specific coexpression networks were constructed from transcriptomic data under different temperature and light conditions. After pruning network edges based on a non-condition-specific literature network and computing and comparing global statistics (such as modularity), a node-wise correlation of network measures (including degree) revealed that similar regulatory pathways were activated under both high-light and cold treatments (Garcia-Molina et al. 2020).

Figure 6:

Comparative analysis of condition-specific stress-response networks. (i) Two networks inferred under distinct stress conditions are shown side-by-side. (ii) Direct comparison applies global and local network statistics (e.g. degree distribution, modularity, centrality) to quantify similarities and differences symbolized by a pair of bar-chart. In parallel, (iii) an alignment step zooms in on three representative nodes from each network: matched node pairs are enclosed by grey ovals and conserved edges are retained, illustrating the mapping between the graphs. (iv) The networks are juxtaposed again, with conserved substructures highlighted and rewired interactions emphasized, allowing direct visualization of shared modules as well as condition-specific changes.

In the context cancer studies this approach has also proven to be effective revealing condition-specific pathway alterations and yielding more biologically interpretable results than simpler expression-based methods. Here, each patient’s gene expression values for metabolic enzymes are used to build a weighted network of metabolites, where edges represent reactions shared by those metabolites. These large networks are then locally aggregated using spectral clustering to form groups of densely connected metabolites (‘supernodes’), thereby reducing complexity and emphasizing pathway-level structure. Next, a specialized Jensen–Shannon divergence-based distance is computed between the networks of different patients. Clustering patients using this network distance matrix often outperforms direct clustering of raw gene expression in distinguishing cancer subtypes, reflecting the additional biological context captured by modeling metabolic interactions (Manipur et al. 2020).

Other methods, such as DINGO and DDN3.0, integrate the condition-specific networks into new networks highlighting rewired or differentially active connections between conditions (Ha et al. 2015; Fu et al. 2024). Both approaches rely on comparing condition-specific dependencies in a joint way, thus identifying a shared “core” along with edges unique to each condition and then merging these differences into a single, newly constructed differential network. This not only pinpoints where network topology diverges under specific stresses but also provides a focused view of likely drivers in stress responses. DDN3.0 was used to analyze proteomic data from 200 human arterial samples linked to early atherosclerosis. Differential comparisons of coronary and aortic proteomes from normal and atherosclerotic tissues revealed significant changes in mitochondrial protein interaction patterns, indicating divergent mitochondrial dynamics. Of 26 identified rewiring hub proteins, tricarboxylic acid cycle components were notably enriched, along with major alterations in necrosis factor, insulin receptor, PPAR-α/γ signaling and rewiring of LXR/RXR activation and collagen system pathways (Fu et al. 2024; Herrington et al. 2018).

When node sets are not fully overlapping across the networks, as is the case for cross-species comparisons or when different experimental platforms yield partially distinct sets of genes or proteins, comparative analysis as described above require additional measures to be applied. One common approach is orthology mapping, where corresponding nodes are matched according to known evolutionary relationships. Tools like OrthoMCL can be used to identify sets of orthologous genes, enabling construction of comparable “orthologous” networks in each organism (Li et al. 2003). Other approaches, such as those in the GRAAL (Graph Aligner) family, directly utilize the structure of networks to be compared by systematically mapping substructures between networks, aiming to align nodes that are topologically similar (Kuchaiev et al. 2010). While the original GRAAL focuses on topological structure, later iterations, including MI-GRAAL, C-GRAAL, and L-GRAAL, can incorporate multiple sources of node similarity such as sequence data (Kuchaiev and Przulj 2011; Malod-Dognin and Pržulj 2015; Memišević and Pržulj 2012).

Also machine learning based methods are used in network alignment, mainly through a mapping of networks to an embedding space, were individual nodes are represented based on their topological structure and attribute information, or in an end-to-end fashion by GNN-based approaches (Chow et al. 2018; Elhesha et al. 2019; Tang et al. 2025).