arXiv:2605.04930v1 [cs.LG] 6 May 2026
When Does Gene Regulatory Network Inference Break? A Controlled Diagnostic Study of Causal and Correlational Methods on Single-Cell Data
Miguel Fernandez-de-Retana∗ Faculty of Engineering University of Deusto Bilbao, Spain [email protected] Unai Zulaika† Faculty of Engineering University of Deusto Bilbao, Spain [email protected]
Rubén Sánchez-Corcuera† Faculty of Engineering University of Deusto Bilbao, Spain [email protected]
Aritz Bilbao-Jayo† Faculty of Engineering University of Deusto Bilbao, Spain [email protected]
Aitor Almeida† Faculty of Engineering University of Deusto Bilbao, Spain [email protected]
Abstract Despite theoretical advantages, causal methods for Gene Regulatory Network (GRN) inference from single-cell RNA-seq data consistently fail to match or outperform correlation-based baselines in many realistic benchmarks, a persistent puzzle which casts doubt on the value of causality for this task. We argue that existing benchmarks are insufficiently controlled to answer this question because they evaluate on real or semi-real data where multiple pathologies co-occur, confounding failure modes, and obscuring the specific conditions under which different inference methods excel or fail. To address this gap, we introduce a controlled diagnostic framework that isolates seven biologically motivated pathologies (dropout, latent confounders, cell-type mixing, feedback loops, network density, sample size, and pseudotime drift) and measure how six representative methods spanning three inference paradigms degrade as each pathology intensifies. Across 6,120 controlled experiments, we find that causal methods genuinely dominate in clean and structurally favorable regimes, but specific pathologies (notably dropout and latent confounders) selectively neutralize their advantages. We further introduce an errortype decomposition that reveals methods with similar aggregate accuracy commit qualitatively different errors. To probe whether single-pathology effects persist when multiple stressors co-occur, we perform an interaction sweep over the three most impactful pathologies and find that their joint effects are sub-additive, while also exposing density-conditional cross-overs invisible to single-dial analysis. Our findings offer a nuanced understanding of when and why different methods succeed or fail for GRN inference, providing actionable insights for method development and practical guidance for practitioners.3 Keywords: Benchmarks, Causal Inference, Computational Biology and Bioinformatics, Network Analysis, Systems Biology ∗ Corresponding author. † Equal contribution. 3 Code available at: https://github.com/miguelfrndz/GRN-Inference-Diagnosis
Preprint.
1
Introduction
Inferring Gene Regulatory Networks (GRNs) from single-cell transcriptomic data is a central problem in computational biology, with profound implications for understanding cellular function, disease mechanisms, drug response, and therapeutic interventions, among others [Kamimoto et al., 2023, Aibar et al., 2017]. Fundamentally, a GRN is a directed graph G = (V, E) encoding regulatory relationships where nodes (V ) represent genes and edges (E) represent functional interactions, such as transcription factors activating or repressing target genes in specific cellular contexts. More formally, the problem of GRN inference can be framed as a structure learning task, where given expression data for p genes across n cells, the goal is to recover the directed adjacency matrix W ∈ Rp×p of the underlying regulatory graph G, where Wij ̸= 0 indicates a regulatory influence of gene i on gene j and the sign of Wij indicates activation or repression. Historically, a wide variety of methods have been proposed for GRN inference, ranging from simple correlation-based approaches (e.g., Pearson correlation, mutual information) to more sophisticated causal inference methods (e.g., PC, GES, NOTEARS) that explicitly model the data-generating process and attempt to recover directed edges. Theoretically, causal methods should have a fundamental advantage: by modeling the underlying generative process, they can, in principle, distinguish direct regulatory interactions from indirect associations and confounding effects, leading to more accurate and interpretable GRN reconstructions. However, a growing body of empirical evidence reveals a puzzling finding: simple correlation-based methods such as Pearson correlation and GENIE3 [HuynhThu et al., 2010] frequently match or outperform theoretically stronger causal inference methods on real single-cell data. This has been documented systematically by CausalBench [Chevalley et al., 2025], which found that causal methods do not reliably outperform baselines on perturbational benchmarks, and by geneRNIB [Nourisa et al., 2025], whose living benchmark reported that simple correlational models often outperformed more complex approaches. These findings have led some to question whether the additional complexity of causal methods is justified for GRN inference at all. We argue that this conclusion is premature. Existing benchmarks evaluate on real or semi-real data where multiple data pathologies, understood as biological or technical challenges that violate the assumptions of different inference methods, co-occur: dropout artifacts, latent confounders, cell-type heterogeneity, and non-stationarity are all present simultaneously. Because these pathologies are confounded with each other and with the unknown ground-truth, benchmarks are underpowered to attribute failures to specific causes. In other words, observing that method A beats method B on dataset D does not tell us why, because we do not know which pathology (or combination of pathologies) is responsible for the performance gap. This lack of diagnostic resolution means we cannot identify the specific conditions under which different methods excel or break, and without knowing why, we cannot improve methods or advise practitioners. Our Approach. We introduce a controlled diagnostic framework inspired by ablation studies in Machine Learning (ML) and perturbation-based analyses to isolate method failure mechanisms under specific pathologies and quantify their impact on performance. We construct a synthetic simulator building upon the previous work by Dibaeinia and Sinha [2020] that generates single-cell expression data from a known linear structural causal model with seven independently controllable pathology “dials”. Each dial corresponds to a specific biological or technical challenge: dropout, latent confounders, cell-type mixing, feedback loops, network density, sample size, and non-stationarity along a pseudotime trajectory (pseudotime drift). By sweeping one dial at a time while holding all others fixed, we isolate the mechanism of failure for each pathology and quantify its impact on performance. This approach allows us to answer questions such as: Does dropout selectively neutralize the advantage of causal methods? Do latent confounders degrade all methods equally? Are some methods more robust to feedback loops than others? By systematically characterizing the degradation patterns of six representative methods spanning three inference paradigms (correlational, tree-ensemble, and causal), we provide a thorough understanding of when and why methods succeed or fail for GRN inference. We evaluate six representative methods spanning five method classes: correlational (Pearson, mutual information), tree-ensemble (GENIE3), constraint-based causal (PC), score-based causal (GES), and optimization-based causal (NOTEARS); using two complementary metrics: Area Under the Precision–Recall Curve (AUPRC) for characterizing the edge-recovery accuracy of each method as pathology intensifies, and a novel error-type decomposition that classifies each predicted edge as true, reversed, confounded, spurious, or missed. 2
Key Findings. Across 2,100 controlled experiments (7 pathologies × 5 levels × 10 replicates × 6 methods) for both linear and nonlinear SCMs, and a complementary interaction sweep over dropout, latent confounders, and density, we find that no method breaks uniformly: causal methods (NOTEARS, GES) dominate under clean conditions (AUPRC > 0.94) and remain the top performer at the hardest level of 6 of the 7 pathology dimensions, but heavy dropout is the exception, where Pearson correlation outperforms all methods. Dropout is also the most discriminating: MI and GENIE3 collapse (∆AUPRC ≈ −0.7), while Pearson degrades gracefully (∆AUPRC = −0.28). Latent confounders and cell-type mixing are the great equalizers, compressing all methods into narrow bands at comparable rates. Error decomposition shows that similar aggregate AUPRC can hide qualitatively different behaviors (e.g., methods may differ in whether their errors are orientation reversals, confounder-induced false positives, or pure hallucinations). Finally, the interaction sweep shows that joint pathologies degrade performance markedly less than the sum of their individual effects, while also revealing density-conditional cross-overs invisible to single-dial analysis.
2
Related Work
GRN Inference Benchmarks. BEELINE [Pratapa et al., 2020] benchmarked 12 algorithms on synthetic (BoolODE) and curated networks, establishing AUROC and AUPRC as standard metrics for GRN evaluation. The DREAM3-5 challenges [Prill et al., 2010, Marbach et al., 2012] systematically evaluated methods on simulated multi-factorial perturbation data. More recently, CausalBench [Chevalley et al., 2025] introduced a large-scale benchmark using real CRISPR single-cell perturbation data with biologically grounded metrics, and geneRNIB [Nourisa et al., 2025] proposed a living benchmark with eight causal-inference metrics across five diverse datasets, with a public leaderboard for continuous evaluation. In both cases, correlation-based methods were found to perform on par with, or better than, more advanced causal methods. Our work is complementary: rather than adding another benchmark on real data, we go in the opposite direction and ask why methods fail, using controlled simulations where the data-generating process is known and specific pathologies can be isolated. GRN Inference Methods. GENIE3 [Huynh-Thu et al., 2010] and its accelerated variant GRNBoost2 [Moerman et al., 2019] remain dominant baselines, forming the backbone of the SCENIC pipeline [Aibar et al., 2017, Bravo González-Blas et al., 2023]. Recent Deep Learning (DL) approaches include DeepSEM [Shu et al., 2021], which parameterizes the GRN as a structural equation model within a Variational Autoencoder (VAE), and single-cell foundation models such as scGPT [Cui et al., 2024], GeneFormer [Theodoris et al., 2023], and scPRINT [Kalfon et al., 2025]. While these methods have shown promising results, a critical evaluation by Ahlmann-Eltze et al. [2025] found that DL-based perturbation prediction does not yet outperform simple linear baselines on well-calibrated metrics. Causal Structure Learning. Constraint-based methods such as Peter–Clark (PC) [Spirtes et al., 2000] test conditional independences to reconstruct a graph skeleton, while score-based methods such as Greedy Equivalence Search (GES) [Chickering, 2002] search over equivalence classes by optimizing a consistent scoring criterion (e.g., BIC). NOTEARS [Zheng et al., 2018] reformulated DAG learning as continuous optimization with an algebraic acyclicity constraint, spawning a family of differentiable methods including DAGMA [Bello et al., 2022], DAG-GNN [Yu et al., 2019], GraN-DAG [Lachapelle et al., 2020], DCDI [Brouillard et al., 2020], DCD-FG [Lopez et al., 2022], and SDCD [Nazaret et al., 2024]. In this work, we focus on the original NOTEARS formulation as a representative optimization-based method, and GES as a representative score-based method, to provide a clear contrast with correlational baselines; our diagnostic contribution is orthogonal to architectural improvements. Perturbation-Based Evaluation. Genome-scale Perturb-seq [Replogle et al., 2022] has made interventional single-cell data widely available at the scale needed to evaluate directed regulatory hypotheses, and has thus become a natural ground-truth for GRN inference evaluation. This has enabled two related but distinct evaluation paradigms: network-inference benchmarks such as CausalBench, which compare inferred edges against perturbational evidence, and perturbationprediction methods such as GEARS [Roohani et al., 2024] and CellOracle [Kamimoto et al., 2023], which evaluate whether a model can predict transcriptional responses under held-out interventions. However, predictive success under perturbation is not equivalent to causal structure recovery: a model may predict average expression shifts without identifying direct regulators, and recent critiques 3
show that complex perturbation predictors can fail to outperform simple linear baselines under well-calibrated metrics [Ahlmann-Eltze et al., 2025]. Conversely, perturbational benchmarks provide realistic validation but still entangle biological heterogeneity, measurement noise, intervention strength, and target-selection biases, which is precisely the entanglement our controlled framework seeks to disentangle.
3
Diagnostic Framework
3.1
Data Generating Process
Our goal is not to simulate every feature of single-cell biology, but to build a diagnostic test bed in which the ground-truth regulatory graph is known and individual failure mechanisms can be switched on independently. We use a linear additive-noise Structural Causal Model (SCM) [Pearl, 2009] as the primary data-generating process and a common baseline for GRN inference. This choice is deliberately favorable to causal discovery methods such as NOTEARS and GES, whose assumptions are closest to this model; consequently, any degradation they exhibit cannot be attributed to an unfairly mismatched simulator. Lastly, to check that our conclusions are not artifacts of linearity, we repeat the full pathology sweep under a nonlinear SCM in Appendix C, where the qualitative degradation patterns are preserved. Ground-Truth Graph. Given p genes, we sample a directed acyclic graph (DAG) G = (V, E) using the natural node order 0 → 1 → · · · → p − 1 as the topological order. For each target gene j, the number of parents is sampled as mj ∼ Binomial(j, ρ), where ρ controls the expected density. Conditional on mj , parents are sampled without replacement from {0, . . . , j − 1} using normalized selection weights qij as in (1): (i + 1)−1 qij = Pj−1 , −1 ℓ=0 (ℓ + 1)
i<j
(1)
so, in practice, earlier genes are more likely to act as regulators. This is a parent-count model with biased parent selection. Edge signs are sampled uniformly from {−1, +1} and magnitudes from U(0.5, 1.0), giving Wij ∈ [−1, −0.5] ∪ [0.5, 1] for selected edges. Linear SCM. Conditioned on the graph, each cell c ∈ {1, . . . , n} is an independent draw from the (c) SCM, where Xj ∈ R is the simulated expression level of gene j in cell c: (c)
Xj
=
X
(c)
Wij Xi
(c)
+ εj ,
(c) iid
εj ∼ N (0, σ 2 )
(2)
i∈PG (j) (c)
where PG (j) denotes the functional parents of gene j in G and εj is the exogenous noise term, modeling unobserved factors and intrinsic stochasticity. The noise variance σ 2 controls the signal-tonoise ratio (SNR) of the data, with higher σ making inference more challenging. In matrix form, taking X ∈ Rn×p as the data matrix of n cells and p genes, W ∈ Rp×p as the weighted adjacency matrix, and ε ∈ Rn×p as the noise matrix, the linear SCM can be expressed as: X = ε (I − W)−1
(3)
where I is the identity matrix. In the default DAG setting this transformation encodes the cumulative effects of all directed paths; when the feedback pathology is enabled, back-edges are added after the DAG is sampled and rescaled if needed for numerical stability, but evaluation still uses the original acyclic graph as the directed ground-truth. 3.2
Pathology Dimensions
Each experiment activates one pathology dial, while keeping the remaining dials at their (benign) defaults. Unless a dial explicitly changes it, the simulator uses p = 25 genes, n = 800 cells, density ρ = 0.1, noise scale σ = 1.0, no dropout, no latent confounders, no cell-type mixing, no feedback, and no pseudotime drift. We define the seven pathology dimensions as follows: 4
P1. Dropout (δ ∈ [0, 0.8]). After expression is generated, each entry is zeroed with an expressiondependent probability, as in (4): h i (c) (c) (c) (c) (c) X̃j = Xj Bj , Bj ∼ Bernoulli 1 − exp −λ (Xj − Xmin ) (4) where Xmin is the global expression minimum across all cells and genes, and λ > 0 is calibrated via (c) binary search so that the marginal dropout rate E[1 − Bj ] = δ for δ ∈ {0, 0.2, 0.4, 0.6, 0.8}. This models the dropout phenomenon in single-cell RNA-seq, where lowly expressed genes are often not detected, leading to zero-inflated data. P2. Latent Confounders (k ∈ {0, 2, 4, 8, 16}). We add k unobserved Gaussian factors to the exogenous noise terms as in (5): (c)
(c)
(c) εj ← εj + L⊤ j z ,
z(c) ∼ N (0, Ik )
(5)
where L ∈ Rp×k is a sparse loading matrix and has entries sampled from N (0, 1) and then masked independently so that each confounder loads on approximately 30% of genes, modeling the effect of unobserved hidden confounding factors. P3. Cell-Type Mixing (α ∈ {0, 0.1, 0.25, 0.4, 0.5}). We draw a fraction (1 − α) of cells from the primary SCM and a fraction α from a second SCM with an independently sampled graph as in (6): XA X= , |XA | = (1 − α)n, |XB | = αn (6) XB where XA ∼ SCM(W1 ) and XB ∼ SCM(W2 ). The mixture is shuffled before inference, and evaluation is performed, as in all cases, against the primary graph W1 . P4. Feedback Loops (ϕ ∈ {0, 0.1, 0.2, 0.3, 0.5}). For each edge in the base DAG, a reverse edge is added with probability ϕ. Back-edge signs are sampled uniformly from {−1, +1} and magnitudes from U(0.1, 0.3), making feedback weaker than the original DAG edges. If the resulting weighted graph has spectral radius at least 0.9, all weights are rescaled by 0.85/ρspec (W) for numerical stability. P5. Network Density (ρ ∈ {0.05, 0.1, 0.15, 0.2, 0.3}). We vary the density parameter used in the parent-count model, mj ∼ Binomial(j, ρ). Thus, ρ controls the expected sparsity of the base DAG and the average number of parents per gene, which in turn affects the complexity of the inference task. Sparser graphs are generally easier to infer, while denser graphs with more parents per gene create more complex dependencies. P6. Sample Size (n ∈ {200, 400, 800, 1600, 3200}). We vary the number of observed cells n, thus isolating the statistical-power effect from changes in graph complexity. P7. Non-Stationary Pseudotime Drift (τ ∈ {0, 0.2, 0.5, 1.0, 1.5}). We partition cells into 10 equal chunks along pseudotime t ∈ [0, 1] and generate each chunk from a time-varying SCM W(t) as in (7) where τ controls the degree of non-stationarity. W(t) = W · 1 + τ (t − 0.5) (7) 3.3
Inference Methods
We evaluate six representative methods spanning correlational, tree-based, constraint-based, scorebased, and continuous-optimization paradigms. Each returns a score matrix S ∈ Rp×p , where larger Sij indicates stronger evidence for i → j. We use default hyperparameters and standard Python implementations to preserve the diagnostic focus on degradation patterns under controlled pathologies, rather than optimizing performance through hyperparameter tuning, which would confound our ability to attribute failures to specific causes. Pearson correlation computes Sij = |corr(X·i , X·j )|, yielding symmetric association scores. Mutual information (MI) estimates pairwise MI after equal-frequency discretization into 6 bins [Cover and Thomas, 1991], also producing symmetric scores. GENIE3 [Huynh-Thu et al., 2010] regresses each target gene on all others using a Random Forest with 50 trees and feature importances as asymmetric edge scores. Peter–Clark (PC) [Spirtes et al., 2000] performs Fisher z conditionalindependence tests up to conditioning-set size 2, returning the skeleton weighted by marginal 5
correlation. Greedy Equivalence Search (GES) [Chickering, 2002] uses a greedy forward BIC search with at most 3 parents per node and a variance-based proxy topological order. NOTEARS [Zheng et al., 2018] solves the penalized linear least-squares objective in (8) with algebraic acyclicity constraint h(W), then thresholds small coefficients and returns absolute edge weights. min W
s.t. 3.4
1 ∥X − XW∥2F + λ∥W∥1 2n h(W) = tr(eW◦W ) − p = 0
(8)
Evaluation Metrics
All methods are evaluated against the original DAG G and under both the linear and nonlinear SCMs. We use two complementary metrics to characterize performance degradation as pathology intensifies. sym AUPRC. For the undirected AUPRC, we collapse edge direction by taking Sij = max(Sij , Sji ) sym and Aij = 1{Aij + Aji > 0}. Precision–recall is then computed over the upper-triangular entries. This is the primary metric for comparing all methods because it does not penalize symmetric methods for lacking orientation. For the directed AUPRC, we compute precision–recall over all off-diagonal ordered pairs, treating (i → j) and (j → i) as distinct predictions. This metric is computed for every method; symmetric methods are not excluded, but their tied scores in both directions naturally limit directional recovery. P Error-Type Decomposition. Given a threshold at the top-K scored edges, where K = ij Aij is the number of ground-truth directed edges, we classify each predicted edge (i, j) into one of five categories relative to the ground-truth DAG G:
True : Reversed : Confounded : Missed : Spurious :
(i, j) ∈ E(G) (j, i) ∈ E(G) but (i, j) ∈ / E(G) i, j share a common ancestor in G but no direct edge (i, j) ∈ E(G) but (i, j) not in top-K none of the above
(9) (10) (11) (12) (13)
This decomposition reveals whether a method fails by reversing directions, hallucinating confounderinduced associations, or producing entirely spurious predictions, distinctions that aggregate AUPRC can mask. Two methods with the same AUPRC may therefore require different methodological fixes.
4
Experiments and Results
We run 7 pathology sweeps × 5 severity levels × 10 replicates × 6 methods, for 2,100 linear SCM experiments, repeat the grid under a nonlinear SCM in Appendix C, and complement these single-dial sweeps with an interaction analysis over the three most informative pathologies in Section 4.4. The main text emphasizes undirected AUPRC because it compares all methods on edge recovery without penalizing symmetric scores for missing orientation; directed AUPRC and runtime results are reported in the Appendix. 4.1
Degradation Curves
Figure 1 shows that the clean setting is not the hard case: under the default linear SCM, structural methods dominate, with NOTEARS at 0.992 undirected AUPRC and GES at 0.944, compared with 0.811–0.895 for Pearson, MI, and GENIE3. The diagnostic value comes from how this ordering changes as each pathology intensifies, revealing strikingly different degradation profiles across methods and pathologies. Baseline Performance. Under clean conditions (leftmost points in each panel), NOTEARS achieves near-perfect recovery (0.992 AUPRC), followed by GES (0.944), PC (0.923), GENIE3 (0.895), MI (0.813), and Pearson (0.811). Explicitly structural methods with directed output (GES, NOTEARS) dominate, while symmetric correlational methods remain competitive above 0.8 AUPRC. Where Each Method Breaks. Table 1 summarizes the change from the easiest to hardest level of each sweep, and three patterns make the failure modes mechanistically explicit. First, dropout breaks 6
0.25 0.00
1.00
0.0
0.2
0.4 0.6 0.8 Dropout Network Density (Edge Density of the Ground-Truth DAG)
0.75 0.50 0.25 0.00
0.05 0.10 0.15 0.20 0.25 0.30 Network Density
0.75 0.50 0.25 0.00
0
1.00
5 10 Confounders Number of Cells (Sample Size)
15
0.75 0.50 0.25 0.00
1000 2000 Number of Cells
3000
Cell-Type Mix Ratio (Cell-Type Mixing)
1.00
AUPRC (Undirected)
0.50
Confounders (Unobserved Latents Driving Gene Groups) AUPRC (Undirected)
0.75
1.00
AUPRC (Undirected)
AUPRC (Undirected)
1.00
AUPRC (Undirected)
AUPRC (Undirected)
AUPRC (Undirected)
Degradation Curves: AUPRC (Undirected) as Each Pathology Intensifies Dropout (Technical Zero-Inflation, Missingness)
0.75 0.50 0.25 0.00
1.00
0.0
0.1
0.2 0.3 0.4 0.5 Cell-Type Mix Ratio Pseudotime Drift (Non-Stationarity Along Pseudotime)
0.75 0.50 0.25 0.00
0.0
0.1
0.2 0.3 Feedback
0.4
0.5
Method
Pearson MI GENIE3 PC GES NOTEARS
0.75 0.50 0.25 0.00
Feedback (Back-Edge Probability, Cycles)
1.00
0.0
0.5 1.0 Pseudotime Drift
1.5
Figure 1: Undirected AUPRC as each pathology intensifies. Lines show mean ± SEM over 10 seeds. Causal methods (GES, NOTEARS) dominate under clean conditions (leftmost points), but degradation is pathology-specific: MI and GENIE3 are most fragile under dropout, confounding and cell-type mixing narrow method differences, and NOTEARS is most robust to feedback, density, sample size, and pseudotime drift. Table 1: AUPRC (undirected) drop from baseline to hardest pathology level. Bold indicates the most robust method per pathology (smallest |∆|); underline indicates the most fragile (largest |∆|). For sample size, positive ∆ (improvement) is expected; bold marks the largest gain. Pathology
Range
Pearson
MI
GENIE3
PC
GES
NOTEARS
Dropout Confounders Cell-type mix Feedback Density Sample size Pseudotime
0 → 0.8 0 → 16 0 → 0.5 0 → 0.5 0.05 → 0.3 200 → 3200 0 → 1.5
−0.28 −0.21 −0.31 −0.18 −0.41 +0.02 −0.14
−0.69 −0.22 −0.29 −0.19 −0.40 +0.04 −0.09
−0.67 −0.25 −0.37 −0.10 −0.34 +0.07 −0.09
−0.41 −0.28 −0.42 −0.12 −0.36 +0.05 −0.20
−0.43 −0.25 −0.41 −0.09 −0.27 +0.01 −0.05
−0.46 −0.34 −0.42 −0.07 −0.05 +0.02 −0.05
discretization- and tree-based scores most aggressively: MI and GENIE3 collapse by 0.69 and 0.67 AUPRC, while Pearson, PC, GES, and NOTEARS converge into a narrow 0.51–0.53 band at δ = 0.8 (Pearson 0.53, NOTEARS 0.53, GES 0.52, PC 0.51). At this level, Pearson is the nominal best by a margin smaller than its standard error, so the operative finding is that heavy dropout neutralizes structural advantages rather than reversing them. Second, latent confounding and cell-type mixing break every method class at comparable rates: at the hardest levels, all six methods fall into narrow bands (0.59–0.70 under confounding; 0.50–0.57 under mixing), so no method class is more robust than the others to these pathologies. Third, NOTEARS is the only method that does not appreciably break under structural pathologies: it remains at 0.923 under feedback, 0.953 under high density, 0.991 at the largest sample size, and 0.944 under pseudotime drift, all settings where correlational approaches lose 0.09–0.41 AUPRC. 4.2
Error-Type Analysis
Aggregate AUPRC hides qualitatively different failure modes. Figure 2 shows the normalized errortype decomposition at the hardest level of each pathology into true, reversed, confounded, spurious, and missed edges. Three structural patterns stand out. Directional Reliability. The reversal rate, i.e., the fraction of top-K predictions that are backwardoriented, is the sharpest separator between structural and correlational methods. NOTEARS and GES maintain reversal rates below 16% across all seven pathologies, while symmetric methods (Pearson, MI, PC) consistently show reversal rates that match their true-positive rates. The contrast is most striking under feedback (ϕ = 0.5): NOTEARS achieves about 70% true and 3% reversed, 7
Normalized Error-Type Composition at the Most Extreme Level
0.4 0.2
0.4 0.2
rson Pea
MI
IE GEN
3
PC
GES
0.0
RS TEA
NO
Network Density = 0.3
1.0
rson Pea
MI
IE GEN
3
PC
GES
0.2
Number of Cells = 3200
Pea
rson
MI
GEN
IE3
PC
GES
NOT
E
ARS
MI
3
IE GEN
PC
GES
RS TEA
NO
0.4
Pea
rson
MI
IE3
GEN
PC
GES
0.4 0.0
S EAR
son
r Pea
NOT
0.0
son
r Pea
MI
GEN
IE3
PC
GES
S EAR
NOT
Edge Type True Reversed Confounded Spurious Missed
0.6 0.2
0.0
0.4
Pseudotime Drift = 1.5
0.8
0.6
0.6 0.2
rson Pea
1.0
0.2
0.0
0.4
NO
Error Share
Error Share
0.4
0.6
0.0
RS TEA
0.8
0.6
0.8
0.2
1.0
0.8 Error Share
0.6
Feedback = 0.5
1.0
0.8 Error Share
0.6
Cell-Type Mix Ratio = 0.5
1.0
0.8 Error Share
Error Share
0.8
0.0
Confounders = 16
1.0
Error Share
Dropout = 0.8
1.0
MI
GEN
IE3
PC
GES
S EAR
NOT
Figure 2: Normalized error-type decomposition at the hardest level of each pathology. Each bar shows the proportion of selected or missed edges assigned to each category. Structural methods recover more true directed edges in several regimes, but their failures are often confident false positives: NOTEARS and GES have more confounded and spurious errors, while Pearson and MI have more missed and reversed edges. NOTEARS maintains the highest True count across most pathologies.
while Pearson scores about 23% true and 23% reversed, a result of its inability to distinguish causal direction. Hallucination Versus Collapse Under Dropout. At δ = 0.8, NOTEARS and GES fill all K prediction slots but with 57% spurious edges, while MI and GENIE3 produce an abundance of spurious edges and near-empty effective prediction sets. Confounded Inflation Under Density. In the dense regime (ρ = 0.3), GES accumulates about 62% confounded predictions at near-zero missed edges, while NOTEARS reaches 53% true, 30% confounded, and 8% reversals. This divergence reflects differing responses to co-expression: GES includes all high-scoring candidate edges and becomes flooded by indirect associations, whereas NOTEARS’s ℓ1 sparsity penalty preferentially selects a smaller subset of edges that are more likely to be direct, thus maintaining a higher true-positive rate at the cost of some missed edges. 4.3
Pathology-Specific Diagnoses
We now synthesize the degradation curves (Figure 1) and error decomposition (Figure 2) into mechanistic explanations for each pathology. Why MI and GENIE3 Collapse Under Dropout. MI relies on equal-frequency discretization into 6 bins, and at δ = 0.8 most observed entries are forced to zero. This severely distorts the entropy estimate, and thus MI collapses to near-random performance. GENIE3’s random forests become similarly uninformative when ∼80% of feature values are zero, as splits cannot reliably separate signals from noise. Pearson correlation, by contrast, is computed over all observations (including zeros) and empirically preserves more of the pairwise edge-score ranking under MNAR zero-inflation, even as magnitudes are attenuated. Why Latent Confounders Degrade All Methods Uniformly. The causal sufficiency assumption, that all common causes of measured variables are also measured, is violated by construction. This is not an implementation failure but a theoretical limitation: without interventional data or instrumental variables, no observational method can distinguish direct causal effects from confounded associations [Pearl, 2009]. The observed compression of hardest-level AUPRC values is therefore expected: confounding removes the clean separation between structural and correlational methods. Why NOTEARS Dominates Under Density and Feedback. NOTEARS combines a linear structural model with an ℓ1 penalty and a global acyclicity constraint. This gives it a strong inductive bias when the target graph is sparse or moderately dense: at ρ = 0.3, it retains AUPRC 0.953 while the next best method, GES, falls to 0.724. Under feedback, the acyclicity constraint h(W) = 0 forces NOTEARS to approximate the best acyclic subgraph of the true (cyclic) structure, which is often close to the original DAG when back-edge magnitudes are small. 8
Best Method per (Dropout, Confounders) Interaction Regime, Grouped by Density Density = 0.1 Density = 0.2
Dropout
Density = 0.05
Density = 0.3
0.8
Pearson 0.67
NOTEARS 0.53
GES 0.40
Pearson 0.40
Pearson 0.55
GES 0.47
Pearson 0.37
Pearson 0.46
Pearson 0.48
GES 0.44
GES 0.44
Pearson 0.40
GES 0.51
GES 0.51
GES 0.50
GES 0.44
0.6
GES 0.88
GES 0.69
GES 0.54
PC 0.57
GES 0.82
GES 0.62
GES 0.53
GES 0.60
GES 0.70
NOTEARS 0.56
NOTEARS 0.54
NOTEARS 0.53
GES 0.63
GES 0.58
GES 0.57
NOTEARS 0.53
0.3
NOTEARS 0.99
NOTEARS 0.90
NOTEARS 0.68
GES 0.66
GES 0.95
NOTEARS 0.83
NOTEARS 0.68
NOTEARS 0.67
NOTEARS 0.89
NOTEARS 0.78
NOTEARS 0.69
NOTEARS 0.68
NOTEARS 0.76
NOTEARS 0.70
NOTEARS 0.70
NOTEARS 0.67
0
NOTEARS 1.00
NOTEARS 0.90
NOTEARS 0.75
GES 0.72
NOTEARS 0.99
NOTEARS 0.91
NOTEARS 0.76
GES 0.70
NOTEARS 0.98
NOTEARS 0.92
NOTEARS 0.80
NOTEARS 0.72
NOTEARS 0.95
NOTEARS 0.88
NOTEARS 0.83
NOTEARS 0.79
0
2
8
16
0
2
8
16
0
2
8
16
0
2
8
16
GENIE3
PC
Confounders
Pearson
Confounders MI
Confounders GES NOTEARS
Confounders
Figure 3: Best method at each (δ, k)-cell, faceted by density ρ. Cells are colored by winning method and annotated with mean AUPRC over seeds. NOTEARS owns the sparse low-dropout corner; Pearson takes over once dropout is heavy; GES surfaces along moderate-confounder, high-dropout at higher density.
4.4
Pathology Interactions
The single-dial sweeps isolate failure mechanisms but say nothing about how methods behave when multiple pathologies act jointly, the norm rather than the exception in real data. To probe this, we run an interaction sweep over the three most informative dials: dropout δ ∈ {0, 0.3, 0.6, 0.8}, latent confounders k ∈ {0, 2, 8, 16}, and network density ρ ∈ {0.05, 0.1, 0.2, 0.3}, yielding 64 joint regimes per method. Figure 3 shows the resulting regime map: at every grid cell, the best-performing method and its mean AUPRC are reported, faceted by density ρ and colored by method class. The full surfaces and per-cell numerical values are reported in Appendix D. Two findings stand out. First, the regime map is structured but not a clean extrapolation of the single-dial ordering: NOTEARS owns the low-dropout half of the grid across all densities, while the heavy-dropout regime is split between Pearson and GES in a density-conditional way. At low density Pearson wins as the single-dial story predicts, but at ρ = 0.3 GES displaces Pearson across the entire δ ≥ 0.6 row, a cross-over invisible to single-dial analysis. Second, joint pathology effects are markedly sub-additive: the AUPRC drop from the clean corner (δ = 0, k = 0, ρ = 0.05) to the diagonally opposite hardest cell is smaller than the sum of the corresponding single-axis drops for every method. As depicted in Table 4 in the Appendix, the gap is largest for MI (−0.96) and GENIE3 (−0.86), already saturated by dropout alone, and smallest for NOTEARS (−0.29).
5
Discussion
From Rankings to Failure Modes. Our controlled experiments provide a mechanistic resolution of benchmark rankings. Different method classes fail for different reasons: causal methods are superior in clean and structurally favorable regimes, dropout selectively damages discretization- and tree-based scores, and confounding compresses differences across all classes. These co-occurring pathologies can make benchmark leaderboards look contradictory, including the familiar cases where simple correlational baselines appear to beat causal methods. Our factorial interaction sweep further shows that the resulting joint failures are largely a sub-additive composition of single-dial mechanisms, so the single-pathology ordering remains a useful reference point. Practical Guidance. Our results yield concrete recommendations: (i) under heavy dropout, prefer Pearson correlation or impute before applying causal methods; (ii) under suspected confounding, use interventional data; (iii) for sparse networks with adequate sample size, NOTEARS provides the best accuracy-per-compute tradeoff (Figure 10 in Appendix E); (iv) always examine error decomposition, not just AUPRC, since qualitatively different errors have different downstream consequences. Limitations. Our study has several limitations that should guide interpretation and future work. (Linear Gaussian SCM) Our primary analysis uses the same linear additive-noise model assumed by NOTEARS and GES. However, the pathologies we study are not artifacts specific to linearity, and Appendix C shows that a nonlinear SCM produces qualitatively similar degradation patterns. (Scale) We use p = 25 genes, representative of a transcription-factor module but far from genome scale, where complexity may change method rankings. (Simulation-first) Our evaluation is entirely synthetic. A systematic study linking simulation pathologies to measurable properties of real perturbational datasets would strengthen the bridge between controlled and real-data evaluation. 9
6
Conclusion
We have presented a controlled diagnostic framework for understanding when and why GRN inference methods fail. By isolating seven biologically motivated pathologies and measuring degradation curves for six representative methods, we provide the first systematic attribution of method failure to specific data properties. Our error-type decomposition reveals that aggregate accuracy metrics mask qualitatively different failure modes with distinct biological implications: structural methods fail by hallucinating high-confidence spurious edges, while correlational and tree-based methods are more prone to collapse, reversal, or under-prediction in other regimes. We find that benchmark rankings are best understood through pathology-specific failure modes rather than overall method labels: the widespread observation that “correlation beats causal” on real data is one consequence of this broader pattern, driven in particular by dropout and latent confounders that selectively neutralize structural advantages. An interaction sweep over the three most informative pathologies further shows that joint effects compose sub-additively and that cross-over effects can arise (e.g., GES displacing Pearson at high density and heavy dropout) that single-dial analysis alone cannot anticipate. Our framework provides a roadmap for future method development and evaluation, allowing researchers to diagnose new methods against specific pathologies, identify pathologies most in need of methodological innovation (latent confounders remain essentially unsolved from observational data alone), and make informed method choices based on their research goals.
Acknowledgments and Disclosure of Funding This work was supported by the Basque Government under grant DEUSTEK5 – Human-Centric Computing for Smart Sustainable Communities and Environments (IT1582-22), and by the European Union’s Horizon Europe research and innovation programme under the PROTECT-CHILD (Grant Agreement No. 101137423) and LATE-AYA (Grant Agreement No. 101214326) projects.
References Constantin Ahlmann-Eltze, Wolfgang Huber, and Simon Anders. Deep-learning-based gene perturbation effect prediction does not yet outperform simple linear baselines. Nature Methods, 22(8): 1657–1661, 2025. Sara Aibar, Carmen Bravo González-Blas, Thomas Moerman, Vân Anh Huynh-Thu, Hana Imrichova, Gert Hulselmans, Florian Rambow, Jean-Christophe Marine, Pierre Geurts, Jan Aerts, et al. SCENIC: Single-cell regulatory network inference and clustering. Nature Methods, 14(11): 1083–1086, 2017. Kevin Bello, Bryon Aragam, and Pradeep Ravikumar. DAGMA: Learning DAGs via M-matrices and a log-determinant acyclicity characterization. In Advances in Neural Information Processing Systems, volume 35, 2022. Carmen Bravo González-Blas, Seppe De Winter, Gert Hulselmans, Nikolai Hecker, Irina Matetovici, Valerie Christiaens, Suresh Poovathingal, Jasper Wouters, Sara Aibar, and Stein Aerts. SCENIC+: Single-cell multiomic inference of enhancers and gene regulatory networks. Nature Methods, 20 (9):1355–1367, 2023. Philippe Brouillard, Sébastien Lachapelle, Alexandre Lacoste, Simon Lacoste-Julien, and Alexandre Drouin. Differentiable causal discovery from interventional data. In Advances in Neural Information Processing Systems, volume 33, pages 21865–21877, 2020. Mathieu Chevalley, Yusuf H Roohani, Arash Mehrjou, Jure Leskovec, and Patrick Schwab. CausalBench: A large-scale benchmark for network inference from single-cell perturbation data. Communications Biology, 8(1):412, 2025. David Maxwell Chickering. Optimal structure identification with greedy search. Journal of Machine Learning Research, 3:507–554, 2002. Thomas M Cover and Joy A Thomas. Elements of Information Theory. Wiley-Interscience, 1991. 10
Haotian Cui, Chloe Wang, Hassaan Maan, Kuan Pang, Fengning Luo, Nan Duan, and Bo Wang. scGPT: Toward building a foundation model for single-cell multi-omics using generative AI. Nature Methods, 21(8):1470–1480, 2024. Payam Dibaeinia and Saurabh Sinha. SERGIO: A single-cell expression simulator guided by gene regulatory networks. Cell Systems, 11(3):252–271, 2020. Vân Anh Huynh-Thu, Alexandre Irrthum, Louis Wehenkel, and Pierre Geurts. Inferring regulatory networks from expression data using tree-based methods. PLoS ONE, 5(9):e12776, 2010. Jérémie Kalfon, Jules Samaran, Gabriel Peyré, and Laura Cantini. scPRINT: Pre-training on 50 million cells allows robust gene network predictions. Nature Communications, 16(1):3607, 2025. Kenji Kamimoto, Blerta Stringa, Christy M Hoffmann, Kunal Jindal, Lilianna Solnica-Krezel, and Samantha A Morris. Dissecting cell identity via network inference and in silico gene perturbation. Nature, 614(7949):742–751, 2023. Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien. Gradient-based neural DAG learning. In International Conference on Learning Representations, 2020. Romain Lopez, Jan-Christian Hütter, Jonathan K Pritchard, and Aviv Regev. Large-scale differentiable causal discovery of factor graphs. In Advances in Neural Information Processing Systems, volume 35, pages 19290–19303, 2022. Daniel Marbach, James C Costello, Robert Küffner, Nicole M Vega, Robert J Prill, Diogo M Camacho, Kyle R Allison, The DREAM5 Consortium, Manolis Kellis, James J Collins, and Gustavo Stolovitzky. Wisdom of crowds for robust gene network inference. Nature Methods, 9(8): 796–804, 2012. Thomas Moerman, Sara Aibar Santos, Carmen Bravo González-Blas, Jaak Simm, Yves Moreau, Jan Aerts, and Stein Aerts. GRNBoost2 and Arboreto: Efficient and scalable inference of gene regulatory networks. Bioinformatics, 35(12):2159–2161, 2019. Achille Nazaret, Justin Hong, Elham Azizi, and David Blei. Stable differentiable causal discovery. In Proceedings of the 41st International Conference on Machine Learning, pages 37413–37445, 2024. Jalil Nourisa, Antoine Passemiers, Marco Stock, Berit Zeller-Plumhoff, Robrecht Cannoodt, Christian Arnold, Alexander Tong, Jason Hartford, Antonio Scialdone, Yves Moreau, Yang Li, and Malte D Luecken. geneRNIB: A living benchmark for gene regulatory network inference. bioRxiv, 2025. doi: 10.1101/2025.02.25.640181. Judea Pearl. Causality: Models, Reasoning, and Inference. Cambridge University Press, 2nd edition, 2009. Aditya Pratapa, Amogh P Jalihal, Jeffrey N Law, Arun Bharadwaj, and T M Murali. Benchmarking algorithms for gene regulatory network inference from single-cell transcriptomic data. Nature Methods, 17(2):147–154, 2020. Robert J Prill, Daniel Marbach, Julio Saez-Rodriguez, Peter K Sorger, Leonidas G Alexopoulos, Xiaowei Xue, Neil D Clarke, Gregoire Altan-Bonnet, and Gustavo Stolovitzky. Towards a rigorous assessment of systems biology models: The dream3 challenges. PLoS ONE, 5(2):e9202, 2010. Joseph M Replogle, Reuben A Saunders, Angela N Pogson, Jeffrey A Hussmann, Alexander Lenail, Alina Guna, Lauren Mascibroda, Eric J Wagner, Karen Adelman, Gila Lithwick-Yanai, et al. Mapping information-rich genotype–phenotype landscapes with genome-scale Perturb-seq. Cell, 185(14):2559–2575, 2022. Yusuf Roohani, Kexin Huang, and Jure Leskovec. Predicting transcriptional outcomes of novel multigene perturbations with GEARS. Nature Biotechnology, 42(6):927–935, 2024. Hantao Shu, Jingtian Zhou, Qiuyu Lian, Han Li, Dan Zhao, Jianyang Zeng, and Jianzhu Ma. Modeling gene regulatory networks using neural network architectures. Nature Computational Science, 1(7): 491–501, 2021. 11
Peter Spirtes, Clark N Glymour, and Richard Scheines. Causation, Prediction, and Search. MIT press, 2000. Christina V Theodoris, Ling Xiao, Anant Chopra, Mark D Chaffin, Zeina R Al Sayed, Matthew C Hill, Helene Mantineo, Elizabeth M Brydon, Zexian Zeng, X Shirley Liu, et al. Transfer learning enables predictions in network biology. Nature, 618(7965):616–624, 2023. Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: DAG structure learning with graph neural networks. In International Conference on Machine Learning, pages 7154–7163. PMLR, 2019. Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. DAGs with NO TEARS: Continuous optimization for structure learning. In Advances in Neural Information Processing Systems, volume 31, 2018.
12
A
Inference Method Details
We provide complete descriptions of the six inference methods evaluated in this study, including their algorithmic details, hyperparameter choices, and key assumptions. Pearson Correlation. Pearson correlation computes the normalized covariance between every pair of genes: Sij = |corr(X·i , X·j )|. The absolute value is taken so that both positive and negative associations contribute equally to the score. The result is a symmetric score matrix with entries in [0, 1]. This method makes no structural assumptions and is relatively robust in our dropout sweep because zero-inflation attenuates many pairwise magnitudes without destroying the entire score ranking. Its main limitation is the inability to distinguish direct edges from indirect associations or to recover edge orientation. Mutual Information (MI). Mutual information is estimated by discretizing each gene’s expression into equal-frequency bins and computing the empirical joint entropy: Sij = H(X·i ) + H(X·j ) − H(X·i , X·j ). We use 6 bins per gene, as in standard GRN benchmarks [Pratapa et al., 2020]. The resulting score matrix is symmetric. MI can in principle detect nonlinear dependencies that Pearson misses, but the discretization step is highly sensitive to zero-inflation: at high dropout, large blocks of tied zero values make the empirical ranks and joint entropies unstable, explaining MI’s catastrophic collapse under the dropout pathology. GENIE3. GENIE3 [Huynh-Thu et al., 2010] casts GRN inference as a feature-importance problem. For each target gene j, a Random Forest regressor is trained on all other genes {X·i }i̸=j to predict X·j , and the feature importances are used as asymmetric edge scores: Sij = importance(X·i → X·j ). We use 50 trees per forest with default Scikit-Learn hyperparameters. The resulting score matrix is asymmetric, with Sij encoding how predictive gene i is for gene j. GENIE3 is robust to nonlinear relationships but, like MI, relies on tree-based splits that become unreliable under heavy dropout. Under dense networks, the predictive signal is spread across many regulators, reducing individual importances. PC (Peter–Clark). The PC algorithm [Spirtes et al., 2000] performs constraint-based causal discovery by iteratively removing edges from the complete graph based on conditional independence tests. We use Fisher z-tests at significance level α = 0.05 with a maximum conditioning-set size of 2. The output is a completed partially directed acyclic graph (CPDAG), from which we extract the skeleton (undirected edges) weighted by marginal Pearson correlation magnitude as the score matrix. This makes the effective score matrix symmetric. PC assumes causal sufficiency (no latent confounders) and faithfulness; violation of either leads to spurious or missing edges. With a conditioning-set limit of 2, PC is computationally efficient but may fail to remove edges in high-density settings where longer adjustment paths are needed. GES (Greedy Equivalence Search). GES [Chickering, 2002] is represented by a compact GES-like forward search using a local linear-Gaussian BIC score. We greedily add parents that improve the target’s BIC, limit the maximum number of parents per node to 3, and restrict candidates using a variance-based proxy topological ordering to enforce acyclicity. The resulting asymmetric score matrix records each selected parent’s positive delta-BIC contribution. This implementation captures the score-based, sparsity-constrained behavior we want to test, but it is not a full two-phase CPDAG search. In dense settings, the BIC score may not strongly discriminate among indirect dependencies, causing confounded associations to be included. NOTEARS. NOTEARS [Zheng et al., 2018] reformulates DAG learning as a continuous optimization problem. It minimizes the penalized least-squares objective: min W
1 ∥X − XW∥2F + λ∥W∥1 2n
s.t.
h(W) = tr(eW◦W ) − p = 0
(14)
where h(W) = 0 is a differentiable algebraic constraint that is zero if and only if W is a DAG. We use λ = 0.05 and threshold learned coefficients at |Wij | < 0.1 to obtain the final edge set. The resulting weight matrix is asymmetric and encodes both direction and magnitude. The ℓ1 penalty promotes sparsity and suppresses indirect associations, providing the strongest inductive bias for direct causal edge recovery in sparse, low-noise networks. The acyclicity constraint causes NOTEARS to prefer the best acyclic approximation when feedback loops are present, which often recovers the dominant causal direction. 13
Correlation vs. Causal Method Families Across Pathologies
0.2 0.0
0.2
0.4 0.6 Dropout Network Density
0.6 0.4 0.2 0.0
0.05 0.10 0.15 0.20 0.25 0.30 Network Density
0.6 0.4 0.2 0.0
0.8
0.8
0.8
1.0
0
5 10 Confounders Number of Cells
15
0.6 0.4 0.2 0.0
0.1
1.0
0.8 0.6 0.4 0.2 0.0
0.8
0.0
1000 2000 Number of Cells
0.5
0.8
0.8 0.6 0.4 0.2 0.0
0.0
0.1
0.2 0.3 Feedback
0.4
0.5
Method Class
0.6
Correlation Tree Causal
0.4 0.2 0.0
3000
0.2 0.3 0.4 Cell-Type Mix Ratio Pseudotime Drift
Feedback
1.0 AUPRC (Undirected)
0.4
Cell-Type Mix Ratio
1.0 AUPRC (Undirected)
0.6
1.0 AUPRC (Undirected)
AUPRC (Undirected)
0.8
0.0
Confounders
1.0
AUPRC (Undirected)
AUPRC (Undirected)
1.0
AUPRC (Undirected)
Dropout
0.0
0.5 1.0 Pseudotime Drift
1.5
Figure 4: Method family comparison: mean ± SEM across methods within each family. The causal family has the highest average AUPRC across all sweeps, with the clearest advantage under density, feedback, and sample-size variation. Table 2: Baseline AUPRC under default settings (p = 25 genes, n = 800 cells, ρ = 0.1, no pathologies), linear Gaussian SCM. Mean over 10 seeds.
AUPRC (undirected) AUPRC (directed)
B
Pearson
MI
GENIE3
PC
GES
NOTEARS
0.811 0.412
0.813 0.414
0.895 0.479
0.923 0.468
0.944 0.911
0.992 0.989
Additional Linear SCM Results
We provide additional analyses of the linear SCM results, including method family comparisons, directed AUPRC degradation curves, and baseline AUPRC values under clean conditions. Method Family Comparison Figure 4 aggregates individual methods into families: correlation (Pearson, MI), tree-ensemble (GENIE3), and causal (PC, GES, NOTEARS). The causal family has the highest mean undirected AUPRC across all seven sweeps, but the margin is pathology-dependent: it is largest under density, feedback, and sample-size variation, and smallest when confounding or cell-type mixing compresses all methods. Notably, the causal family’s advantage under the density sweep is large: at ρ = 0.3, the causal family mean is 0.77 AUPRC while the correlation family drops to ≈ 0.55. Conversely, at peak dropout (δ = 0.8), MI’s collapse pulls the correlation mean below 0.4, obscuring the fact that Pearson alone remains competitive. Disaggregating by method (Figure 1) is therefore essential for interpreting family-level aggregates. Directed AUPRC Degradation Curves All main-text analyses use undirected AUPRC to compare edge recovery without penalizing symmetric score matrices. Figure 5 repeats the degradation analysis with directed AUPRC. The main difference is expected: Pearson, MI, and PC lose roughly half their clean AUPRC because their scores are symmetric (producing tied scores in both directions), while GES and NOTEARS retain directed AUPRC close to their undirected values under clean conditions (0.911 and 0.989, respectively; see Table 2). Importantly, the relative degradation patterns across pathologies are preserved: dropout still collapses MI and GENIE3, confounding compresses method differences, and NOTEARS dominates under the density and feedback pathologies. The absolute gap between structural and correlational methods is wider in the directed setting, as edge orientation is precisely what symmetric methods cannot recover. Figure 6 shows the corresponding family-level directed comparison, where the causal family retains a larger margin over the correlation family due to GES and NOTEARS’s asymmetric parent-to-target scoring. Baseline AUPRC Values Table 2 reports the AUPRC values under clean conditions (default pathology levels) for reference. These baseline values establish the upper bound on performance for each method and represent the regime most favorable to causal discovery: a linear Gaussian SCM with adequate sample size, no dropout, no latent confounders, and moderate network density. 14
Degradation Curves: AUPRC (Directed) as Each Pathology Intensifies
0.0
0.4 0.6 0.8 Dropout Network Density (Edge Density of the Ground-Truth DAG)
0.75 0.50 0.25 0.00
0.50 0.25 0
5 10 Confounders Number of Cells (Sample Size)
1.00
0.05 0.10 0.15 0.20 0.25 0.30 Network Density
0.50 0.25 1000 2000 Number of Cells
0.25
1.00
0.75
0.00
0.50
0.1
0.2 0.3 0.4 0.5 Cell-Type Mix Ratio Pseudotime Drift (Non-Stationarity Along Pseudotime)
0.75 0.50 0.25 0.00
0.0
0.1
0.2 0.3 Feedback
0.4
0.5
Method
Pearson MI GENIE3 PC GES NOTEARS
0.75 0.50 0.25 0.00
3000
0.0
Feedback (Back-Edge Probability, Cycles)
1.00
0.75
0.00
15
Cell-Type Mix Ratio (Cell-Type Mixing)
1.00
AUPRC (Directed)
0.75
0.00
0.2
AUPRC (Directed)
0.25
Confounders (Unobserved Latents Driving Gene Groups)
AUPRC (Directed)
0.50
1.00 AUPRC (Directed)
AUPRC (Directed)
0.75
0.00
1.00
AUPRC (Directed)
AUPRC (Directed)
1.00
Dropout (Technical Zero-Inflation, Missingness)
0.0
0.5 1.0 Pseudotime Drift
1.5
Figure 5: Directed AUPRC degradation across all seven pathologies. Symmetric score matrices (Pearson, MI, and the PC skeleton implementation) lose orientation information and operate near their random baseline in directed evaluation, whereas GES and NOTEARS maintain strong directed performance under clean and moderately perturbed conditions. Correlation vs. Causal Method Families Across Pathologies
0.2 0.0
0.2
0.4 0.6 Dropout Network Density
0.6 0.4 0.2 0.0
0.05 0.10 0.15 0.20 0.25 0.30 Network Density
0.6 0.4 0.2 0.0
0.8
0.8
0.8
AUPRC (Directed)
0.4
1.0
0
5 10 Confounders Number of Cells
0.4 0.2
1.0
0.6 0.4 0.2 1000 2000 Number of Cells
0.6
0.0
0.2 0.3 0.4 Cell-Type Mix Ratio Pseudotime Drift
0.5
0.8
0.8 0.6 0.4 0.2 0.0
0.0
0.1
0.2 0.3 Feedback
0.4
0.5
Method Class
0.6
Correlation Tree Causal
0.4 0.2 0.0
3000
0.1
Feedback
1.0
0.8
0.0
15
0.8
0.0
Cell-Type Mix Ratio
1.0
AUPRC (Directed)
0.6
1.0 AUPRC (Directed)
AUPRC (Directed)
0.8
0.0
Confounders
1.0
AUPRC (Directed)
AUPRC (Directed)
1.0
AUPRC (Directed)
Dropout
0.0
0.5 1.0 Pseudotime Drift
1.5
Figure 6: Directed AUPRC method family comparison. The causal family advantage is substantially larger in directed evaluation because GES and NOTEARS produce asymmetric parent-to-target scores, while correlation and tree-ensemble methods either produce symmetric scores or inherit directional limitations from their feature-importance formulation.
The large gap between undirected and directed AUPRC for Pearson, MI, and PC reflects their symmetric score matrices: without orientation signal, these methods are limited to approximately chance-level directed AUPRC relative to the number of ordered pairs. GENIE3’s moderate directed AUPRC (0.479) reflects partial asymmetry from its feature-importance formulation, but this is far below its undirected performance because importance scores do not reliably encode causal direction. GES and NOTEARS maintain near-perfect directed AUPRC (> 0.9), demonstrating that their asymmetric parent-to-target scores recover edge orientation reliably under clean conditions. This clean-condition baseline is the reference against which all degradation in the pathology sweeps should be interpreted.
C
Nonlinear SCM Results
To test whether the diagnostic conclusions depend on the linear Gaussian SCM, we replicate the full sweep under a nonlinear SCM where each gene’s expression is a nonlinear function of its parents plus noise. We choose a simple nonlinearity, the hyperbolic tangent (tanh), which preserves the additive noise structure while introducing bounded nonlinear interactions that are common in gene regulation. 15
Linear vs. Nonlinear SCM: AUPRC (Undirected) Degradation Curves Dropout
AUPRC (Undirected)
1.0
Cell-Type Mix Ratio
Feedback
0.8 0.6 0.4 0.2 0.0
0.0
0.2
0.4
0.6
0.8
0
2
Network Density
1.0 AUPRC (Undirected)
Confounders
4
8
16
0.0
0.1
Number of Cells
0.25
0.4
0.5
0.1
0.2
0.3
0.5
Method / SCM Type
Pearson MI GENIE3 PC GES NOTEARS Linear SCM Nonlinear SCM (Tanh)
0.8 0.6 0.4 0.2 0.0
0.0
Pseudotime Drift
0.05
0.1
0.15
0.2
0.3
200
400
800
1600
3200
0.0
0.2
0.5
1.0
1.5
Figure 7: Linear (solid) vs. nonlinear tanh (dashed) SCM, using undirected AUPRC. The nonlinear SCM changes absolute difficulty and some method rankings, but the diagnostic stressors remain visible: dropout still collapses MI and GENIE3, confounding remains broadly damaging, and NOTEARS remains strongest under high density. Definition 1 (Nonlinear SCM (Tanh)). The nonlinear variant replaces the linear structural assignment with: X (c) (c) (c) (c) iid Xj = tanh Wij Xi + εj , εj ∼ N (0, σ 2 ) (15) i∈PG (j)
For acyclic graphs, sampling proceeds by ancestral ordering. For cyclic graphs (with feedback), we use fixed-point iteration X(t+1) = tanh(X(t) W)+ε, which converges when ∥W∥2 < 1 (guaranteed by our rescaling) since the tanh nonlinearity is 1-Lipschitz. Figure 7 overlays linear (solid) and nonlinear (dashed) degradation curves. The nonlinear sweep should not be read as an identical replay of the linear results. Because the tanh nonlinearity bounds expression values, several settings become easier in absolute undirected AUPRC (e.g., NOTEARS reaches near-perfect AUPRC under clean conditions in both regimes), while strong latent confounding becomes harder for all methods due to the increased expressiveness of nonlinear confounding pathways. Critically, the qualitative ordering of pathology severity is preserved: dropout remains the most destructive pathology for MI and GENIE3 (whose discretization-based scores are equally disrupted under either SCM); confounding compresses method differences regardless of the functional form; and NOTEARS retains the largest advantage in the high-density regime, where its sparsityinducing ℓ1 penalty and acyclicity constraint continue to suppress indirect associations. The nonlinear results thus corroborate the central diagnostic conclusions and demonstrate that our findings are not artifacts of the linear Gaussian modeling assumption. Figure 8 shows the corresponding directed comparison, where the gap between structural and symmetric methods is further amplified because the nonlinearity does not change the fundamental inability of Pearson, MI, and PC to recover edge orientation. Table 3 reports the corresponding baseline AUPRC under the nonlinear SCM at the same default settings. Compared to the linear baseline (Table 2), the bounded tanh nonlinearity increases clean undirected AUPRC for several methods, while NOTEARS retains near-identical performance. Directed AUPRC patterns are preserved: GES and NOTEARS maintain a large gap over symmetric methods, with GENIE3’s directed score improving due to the saturating nonlinearity making parent contributions more distinguishable from indirect associations.
D
Pathology Interaction Sweeps
This appendix expands on the factorial interaction analysis introduced in Section 4.4. We jointly vary the three most informative single-dial pathologies, dropout δ ∈ {0, 0.3, 0.6, 0.8}, latent confounders k ∈ {0, 2, 8, 16}, and network density ρ ∈ {0.05, 0.1, 0.2, 0.3}, on a fully crossed grid (64 joint 16
Linear vs. Nonlinear SCM: AUPRC (Directed) Degradation Curves Dropout
AUPRC (Directed)
1.0
Confounders
Cell-Type Mix Ratio
Feedback
0.8 0.6 0.4 0.2 0.0
0.0
0.2
1.0 AUPRC (Directed)
0.4
0.6
0.8
0
2
Network Density
4
8
16
0.0
0.1
Number of Cells
0.25
0.4
0.5
0.1
0.2
0.3
0.5
Method / SCM Type
Pearson MI GENIE3 PC GES NOTEARS Linear SCM Nonlinear SCM (Tanh)
0.8 0.6 0.4 0.2 0.0
0.0
Pseudotime Drift
0.05
0.1
0.15
0.2
0.3
200
400
800
1600
3200
0.0
0.2
0.5
1.0
1.5
Figure 8: Linear vs. nonlinear SCM comparison using directed AUPRC. GES and NOTEARS retain the strongest orientation signal in clean and moderately perturbed settings, while symmetric methods remain limited by tied scores in both directions. The directed gap between structural and correlational methods is wider than in the undirected case, as orientation recovery is the primary differentiator between method classes. Table 3: Baseline AUPRC under default settings (p = 25 genes, n = 800 cells, ρ = 0.1, no pathologies), nonlinear tanh SCM. Mean over 10 seeds.
AUPRC (undirected) AUPRC (directed)
Pearson
MI
GENIE3
PC
GES
NOTEARS
0.935 0.479
0.925 0.470
0.967 0.591
0.979 0.501
0.960 0.891
0.991 0.983
regimes), with all remaining pathology dials held at their benign defaults and 5 random seeds per cell. The fixed simulator settings match the single-dial sweeps so that interaction-specific effects are not confounded with scale. Per-Method Failure Surfaces. Figure 9 shows, for each method, the mean undirected AUPRC over the dropout × confounders plane at each density slice. Three patterns are visible. First, the dropout axis is the dominant degrader for MI and GENIE3 across every density; their surfaces collapse uniformly along this axis with little additional movement once δ > 0.6. Second, NOTEARS retains a high-AUPRC region at low dropout across confounder levels, but this region contracts once dropout reaches δ = 0.6 and is largely lost at δ = 0.8. Third, density modulates absolute AUPRC for every method except NOTEARS, whose ridge stays near 1.0 across the entire density axis; PC and GES also respond to density less than the correlation- and tree-based methods, consistent with Table 1. Sub-Additivity of Joint Pathologies. Table 4 quantifies the sub-additivity claim from Section 4.4. For each method we compute the AUPRC drop from the cleanest joint cell (δ = 0, k = 0, ρ = 0.05) to the worst joint cell (δ = 0.8, k = 16, ρ = 0.3), and compare it against the sum of the corresponding single-dial drops. The interaction term, defined as ∆joint − ∆additive , is negative for every method, indicating that pathologies share failure modes rather than compounding independently. The effect is largest for MI and GENIE3 (already saturated by dropout alone) and smallest for NOTEARS, whose degradation is closest to additive. Reading the Winner Map. The regime map in Figure 3 (Section 4.4) summarizes the surfaces by reporting only the best method at each cell. NOTEARS owns the low-dropout half of the grid across all densities, with the largest margin in the sparser slices. The heavy-dropout regime is more nuanced: at low density Pearson wins, consistent with its graceful degradation under MNAR zero-inflation, but as density increases GES progressively reclaims those cells, and at ρ = 0.3 GES displaces Pearson entirely across the δ ≥ 0.6 rows. This cross-over is the most visible deviation from the single-dial ordering and is consistent with GES’s forward-search BIC selection benefiting from richer dependency structure once edges are abundant. 17
Joint Pathology Failure Surfaces (Dropout x Confounders, Grouped by Density)
NOTEARS Dropout
GES Dropout
Density = 0.1
Density = 0.2
Density = 0.3
0.8
0.67
0.45
0.36
0.40
0.55
0.46
0.37
0.46
0.48
0.42
0.41
0.40
0.47
0.47
0.47
0.42
0.6
0.87
0.65
0.47
0.57
0.76
0.56
0.51
0.57
0.60
0.48
0.48
0.46
0.54
0.50
0.50
0.47
0.3
0.94
0.71
0.56
0.61
0.81
0.61
0.58
0.59
0.63
0.51
0.52
0.50
0.55
0.52
0.53
0.51
0
0.95
0.70
0.55
0.63
0.82
0.62
0.58
0.59
0.63
0.51
0.52
0.49
0.55
0.51
0.52
0.50
0.8
0.07
0.07
0.06
0.07
0.12
0.12
0.12
0.14
0.21
0.23
0.23
0.18
0.31
0.31
0.34
0.30
0.6
0.09
0.09
0.09
0.07
0.17
0.15
0.14
0.17
0.21
0.25
0.23
0.27
0.33
0.37
0.36
0.34
0.3
0.88
0.65
0.51
0.58
0.78
0.61
0.58
0.57
0.60
0.50
0.51
0.49
0.55
0.51
0.52
0.50
0
0.95
0.70
0.54
0.62
0.82
0.63
0.59
0.58
0.62
0.51
0.52
0.49
0.55
0.51
0.52
0.49
0.8
0.14
0.16
0.10
0.15
0.24
0.22
0.20
0.20
0.30
0.29
0.29
0.31
0.37
0.38
0.39
0.36
0.6
0.71
0.56
0.38
0.49
0.66
0.54
0.47
0.50
0.61
0.49
0.45
0.45
0.54
0.52
0.51
0.47
0.3
0.94
0.79
0.56
0.64
0.85
0.70
0.60
0.60
0.71
0.59
0.57
0.55
0.61
0.58
0.59
0.56
0
0.97
0.82
0.58
0.66
0.88
0.73
0.62
0.65
0.76
0.62
0.60
0.57
0.64
0.60
0.59
0.58
0.8
0.67
0.44
0.36
0.40
0.54
0.46
0.36
0.46
0.47
0.43
0.40
0.40
0.48
0.46
0.48
0.42
0.6
0.87
0.65
0.47
0.57
0.76
0.56
0.51
0.58
0.60
0.48
0.48
0.46
0.56
0.51
0.51
0.48
0.3
0.94
0.72
0.56
0.61
0.83
0.63
0.60
0.60
0.64
0.54
0.54
0.51
0.59
0.55
0.57
0.53
0
0.97
0.77
0.57
0.66
0.91
0.74
0.61
0.63
0.71
0.59
0.57
0.51
0.63
0.60
0.58
0.55
0.8
0.54
0.43
0.40
0.33
0.55
0.47
0.36
0.44
0.48
0.44
0.44
0.40
0.51
0.51
0.50
0.44
0.6
0.88
0.69
0.54
0.57
0.82
0.62
0.53
0.60
0.70
0.55
0.53
0.50
0.63
0.58
0.57
0.51
0.3
0.98
0.85
0.63
0.66
0.95
0.75
0.67
0.65
0.80
0.66
0.64
0.57
0.71
0.66
0.66
0.58
0
0.99
0.86
0.68
0.72
0.95
0.85
0.72
0.70
0.86
0.73
0.69
0.62
0.73
0.69
0.65
0.63
0.8
0.53
0.53
0.37
0.36
0.51
0.44
0.28
0.41
0.37
0.39
0.37
0.39
0.42
0.42
0.45
0.44
0.6
0.62
0.56
0.48
0.55
0.67
0.60
0.51
0.56
0.67
0.56
0.54
0.53
0.57
0.55
0.55
0.53
0.3
0.99
0.90
0.68
0.63
0.93
0.83
0.68
0.67
0.89
0.78
0.69
0.68
0.76
0.70
0.70
0.67
0
1.00
0.90
0.75
0.67
0.99
0.91
0.76
0.63
0.98
0.92
0.80
0.72
0.95
0.88
0.83
0.79
0
2
8
16
0
2
8
16
0
2
8
16
0
2
8
16
Confounders
Confounders
Confounders
Confounders
1.0
0.8
0.6 AUPRC (Undirected)
PC Dropout
GENIE3 Dropout
MI Dropout
Pearson Dropout
Density = 0.05
0.4
0.2
0.0
Figure 9: Per-method failure surfaces over the dropout × confounders plane, faceted by density. Each cell reports mean undirected AUPRC over 5 seeds. Shared [0, 1] colour scale. The dropout axis dominates degradation for MI and GENIE3, while NOTEARS retains strong performance in the low-dropout regime before dropping sharply at higher dropout.
E
Experimental Details
Reproducibility. All experiments use NumPy’s default random number generator with seeds {0, 1, . . . , 9}. The full experiment suite can be reproduced with a single command and run on a standard CPU with consumer-grade specifications. Code and data will be made publicly available upon publication. 18
Table 4: Joint vs. additive AUPRC degradation across the dropout × confounders × density grid. ∆joint is the AUPRC drop from the cleanest to the worst joint cell; ∆d , ∆k , ∆ρ are the corresponding single-dial drops along each axis (other dials held at their cleanest level); ∆add = ∆d + ∆k + ∆ρ ; interaction = ∆joint − ∆add . All entries are AUPRC differences (lower means more degradation); negative interaction values indicate sub-additive composition. Method
AUPRC0
∆joint
∆d
∆k
∆ρ
Interaction
Pearson MI GENIE3 PC GES NOTEARS
0.954 0.949 0.971 0.972 0.989 1.000
0.532 0.646 0.613 0.553 0.549 0.563
0.282 0.878 0.829 0.304 0.444 0.470
0.326 0.326 0.313 0.309 0.267 0.328
0.402 0.400 0.328 0.340 0.259 0.052
−0.477 −0.958 −0.857 −0.399 −0.421 −0.287
Graph Construction. Ground-truth DAGs are constructed with node ordering 0 < 1 < · · · < p − 1 as the topological order. Full details are provided in Section 3.1. Software. Implementation uses NumPy, SciPy, Scikit-Learn (random forests for GENIE3), NetworkX (graph operations), and Matplotlib (visualization). No GPU is required. Accuracy–Runtime Tradeoff. Figure 10 presents the accuracy–runtime Pareto frontier, averaged across all linear-SCM experiments. Runtimes span more than four orders of magnitude: Pearson correlation completes in microseconds (∼ 10−4 s) while GENIE3’s random-forest fitting takes ∼ 1.5s. NOTEARS occupies a compelling middle ground, achieving the highest average undirected AUPRC (0.912 averaged across all pathology levels) at a moderate cost of ∼ 0.16s per experiment. GES is the most efficient structural method, retaining strong directed AUPRC under clean conditions at roughly one-tenth of NOTEARS’s runtime. The Pareto analysis reinforces the practical recommendation from the main text: for sparse networks with adequate sample size and without heavy dropout or confounding, NOTEARS provides the best accuracy-per-compute tradeoff. For settings where runtime is the primary constraint, GES is a strong alternative. Pearson correlation, despite being the fastest method by a wide margin, is dominated in accuracy by all structural methods under clean conditions, and should only be preferred when dropout is heavy (Section 4.3).
Pareto View of Accuracy vs. Runtime Across All Experiments 1.0
AUPRC
0.8 0.6 0.4 0.2 0.0
10 4 Pearson MI
10 1 10 3 10 2 Mean Runtime per Call (s, Log Scale) GENIE3 PC
GES NOTEARS
100
Undirected Directed
Figure 10: Pareto view of accuracy vs. runtime (log scale), averaged across all linear-SCM experiments. Circles: undirected AUPRC; Squares: directed AUPRC. Pearson is fastest (∼ 10−4 s) but less accurate; GES is both fast (∼ 2 × 10−2 s) and structurally accurate; NOTEARS gives the highest average AUPRC at moderate cost (∼ 1.6 × 10−1 s); GENIE3 is slowest (∼ 1.5s) without a commensurate accuracy gain. 19