1 Introduction↩︎

Gene regulatory networks (GRNs) dictate cellular identity and physiological function by encoding the complex interactions between transcription factors (TFs) and target genes [1]–[4]. While the advent of single-cell RNA sequencing (scRNA-seq) has fundamentally advanced GRN inference by resolving bulk-level signal averaging [5], the resulting data is notoriously heterogeneous and sparse [6], [7]. To tackle these challenges, recent progress has shifted from classical statistical methods to graph neural networks (GNNs), particularly graph autoencoder (GAE) architectures [8]–[13]. Despite these architectural advances, we argue that computational GRN inference suffers from a fundamental misalignment with practical biological discovery in two critical dimensions:

First, there is an evaluation misalignment. Standard benchmarks such as BEELINE [14] evaluate GRN inference under transductive settings, where all genes are observed during training. However, discovering regulatory mechanisms for novel modules or rare cell types inherently requires inductive generalization to unseen genes [15], [16]. Moreover, these benchmarks formalize GRN inference as a binary classification task measured by AUROC and AUPRC. In practice, however, biologists rely on AI models to prioritize a small, high-confidence set of novel interactions for costly experimental validation. Global metrics such as AUROC are easily dominated by low-ranked edges, making them poor proxies for the top-\(K\) ranking quality that ultimately drives scientific discovery.

Second, there exists a gene-level inductive generalization gap. Under realistic inductive settings, existing GAE-based methods rely heavily on predefined prior graphs. When holdout genes become disconnected during inference, they lose critical contextual signals and suffer severe performance degradation (Figure 1a). Recent diffusion-based methods partially alleviate this topological isolation, but still fail to generalize effectively to unseen genes, as the highly heterogeneous and sparse nature of scRNA-seq data prevents reliable feature correspondence and can lead to feature collapse.

More fundamentally, current approaches fail to capture the conditional nature of gene regulations. Regulatory interactions are not isolated events, but are conditioned on the paired states of regulator and target genes. Applying a uniform transition to all edges during diffusion, regardless of the activity states of their endpoint genes, violates this fundamental principle of GRNs.

To address these limitations, we introduce a new evaluation paradigm alongside a novel generative framework. First, we propose BEELINE-KGC, a systematic extension of BEELINE that enforces an Inductive Gene-Holdout Split. By recasting evaluation as a Knowledge Graph Completion (KGC) task [17], BEELINE-KGC uses Hits@K and Mean Reciprocal Rank (MRR) to directly measure the top-\(K\) ranking consistency of regulatory predictions. Under this rigorous protocol, existing state-of-the-art models exhibit severe generalization bottlenecks.

On the modeling side, we propose CoDiffGRN, the first gene-regulation co-evolutionary discrete diffusion model for GRN inference. Built upon the D3PM framework [18], CoDiffGRN departs from prior-dependent paradigms and is designed for robust inductive generalization. To mitigate feature collapse on unseen genes, we introduce a Cell-cluster-based Discretization strategy that transforms scRNA-seq profiles into biologically coherent discrete states, allowing unseen genes to inherit prior distributional structure during diffusion.

To capture the conditional nature of gene regulation, we introduce a joint node-edge diffusion process that directly models the co-evolution of discretized gene expression states and network topology at the level of diffusion dynamics. By devising a gene-state-dependent transition matrix, edge transitions are explicitly conditioned on the states of their endpoint genes, promoting the coherent neighborhoods required for top-\(K\) regulatory discovery. Furthermore, to enable scalable generative training under data scarcity, we introduce TF-ALL Subgraph Sampling (TASS), a data augmentation strategy that preserves TF-target integrity while substantially improving scalability and generalization.

Comprehensive experimental results on BEELINE-KGC show that CoDiffGRN achieves state-of-the-art performance and consistently outperforms existing models. Extensive ablation studies examine the contributions of the proposed gene-regulation co-evolutionary discrete diffusion and TASS framework. These results collectively support the effectiveness of our approach in improving KGC-based ranking performance, enabling more accurate and reliable biological discovery.

2 Related Work↩︎

Evaluation protocols for GRN inference. BEELINE [14] established the standard for GRN inference, primarily utilizing transductive edge-centric splits and global metrics such as AUROC and AUPRC [19]. To improve discriminative sensitivity, negative sampling has evolved from random selection [20], [21] to heuristic strategies [10], [13]. However, these protocols exhibit a fundamental misalignment with practical biological discovery: transductive settings neglect the imperative for generalization to unseen genes, while global metrics fail to prioritize the top-\(K\) candidates required for experimental validation. In this work, we introduce BEELINE-KGC to bridge this gap by recasting evaluation as an inductive KGC task, utilizing ranking-based metrics (Hits@K, MRR) to assess top-\(K\) regulatory predictions.

Gene regulatory network inference methods. Early GRN inference relied on statistical associations like GENIE3 [22] and GRNBoost2 [20], or feature-based deep learning models such as CNNC [21] and DeepSEM [23]. Recently, graph representation learning has become the dominant paradigm, with GAE-based frameworks like GENELink [8], GENELink+ [10] and GNNLink [9] utilizing prior regulatory graphs, while methods like GCLink [12] and GRNFormer [11] introduce contrastive learning and VGAE-based frameworks. More recently, generative models such as RegDiffusion [24], DigNet [25], and Planet [26] have applied continuous and discrete diffusion processes [18], [27] to model the evolution nature of gene regulations.

However, existing methods face critical limitations: traditional approaches struggle with noisy single-cell data, while graph-based models often hinge on prior-graph connectivity. Furthermore, current generative frameworks struggle with feature collapse in inductive scenarios and fail to capture the essential co-evolution between gene and regulation with severe data scarcity. In this work, we propose CoDiffGRN, which employs cell-cluster discretization and joint discrete diffusion to model regulatory dynamics, alongside TASS specifically designed to address the challenges of data scarcity.

3 Rethinking for Practical Discovery: The BEELINE-KGC Benchmark↩︎

3pt

We introduce BEELINE-KGC, a comprehensive benchmark designed to bridge the gap between computational GRN inference and realistic biological discovery. Built upon the BEELINE [14] (CC BY-NC 4.0) collection of experimental scRNA-seq datasets, BEELINE-KGC establishes a rigorous evaluation framework tailored to real-world laboratory constraints.

3.1 Formalization and Problem Setting↩︎

We formalize a gene regulatory network as a directed graph \(\mathcal{G}=(\mathcal{V}, \mathcal{E})\), where nodes \(v \in \mathcal{V}\) represent genes and edges \((u, v) \in \mathcal{E}\) denote regulatory interactions. We partition \(\mathcal{V}\) into TFs (\(\mathcal{V}_{TF}\)) and target genes (\(\mathcal{V}_{tgt}\)). Within the TFs, we further distinguish source TFs (\(\mathcal{V}_{src} \subseteq \mathcal{V}_{TF}\)), defined as regulators with at least one observed outgoing edge in the ground truth.

The corresponding scRNA-seq data is represented by an expression matrix \(\mathbf{X} \in \mathbb{R}^{N \times F_X}\), where \(\mathbf{X}_{i,j}\) denotes the expression level of gene \(v_i \in \mathcal{V}\) in cell \(j \in \{1, \dots, F_X\}\). Given the expression matrix \(\mathbf{X}\) and a subset of observed interactions \(\mathcal{E}_{\text{obs}}\), the primary objective of GRN inference is to recover the unobserved regulatory structure by predicting potential links in \(\mathcal{E} \setminus \mathcal{E}_{\text{obs}}\).

3.2 A Unified Framework for Practical GRN Discovery↩︎

Standard GRN inference benchmarks adopt a transductive setting and rely on global metrics such as AUPRC. However, this paradigm fundamentally diverges from the practical demands of real biological discovery, which requires inductive generalization to novel regulators and high top-\(K\) ranking quality under limited experimental budgets.

Pillar 1: Inductive gene-holdout split. In realistic discovery scenarios, researchers predominantly focus on characterizing novel or unannotated TFs by identifying their downstream regulatory targets, necessitating an inductive setting where models must generalize to unseen regulators [15], [16], [28]–[33]. To simulate this, we extract a subset of \(\mathcal{V}_{src}\) as holdout TFs (\(\mathcal{V}_{\mathrm{hold}}\)), introducing an inductive protocol:

1. Training: The model observes only the subgraph induced by \(\mathcal{V} \setminus \mathcal{V}_{hold}\). Expression profiles and all edges involving holdout TFs are withheld.

2. Evaluation: The model receives the expression profiles for the full gene set \(\mathcal{V}\) and predicts regulations across the entire network. Performance is evaluated solely on edges involving at least one holdout TF.

Pillar 2: KGC reformulation and multi-scale ranking protocol. Practical biological discovery is governed by a discrete experimental budget \(K\), where true interactions are useful only if ranked within the top-\(K\) and lower ranks are effectively equivalent to being undiscovered. However, global metrics like AUPRC are insensitive to top-ranked precision, and are further distorted by the severe class imbalance in GRNs leaving many candidate pairs unevaluated under standard negative sampling.

To resolve this, we reformulate GRN inference as a Knowledge Graph Completion (KGC) problem. Each regulation is modeled as a structural triple \((v_i, \texttt{regulates}, v_j)\), where \(v_i\) is the source TF (head) and \(v_j\) is the target gene (tail). The task of finding targets for a holdout TF naturally translates to tail entity prediction: given a query \((v_i, \texttt{regulates}, ?)\), the model must rank all candidate tails \(t \in \mathcal{V}\) such that the ground-truth targets are placed as high as possible.

Aligned with this KGC formulation, we evaluate GRN inference as a resource-constrained ranking problem, adopting Hits@\(K\) for top-\(K\) recovery under a fixed validation budget and Mean Reciprocal Rank (MRR) for capturing overall ranking quality across heterogeneous budgets: \[\text{Hits@}K = \frac{1}{|\mathcal{Q}|}\sum_{q \in \mathcal{Q}} \mathbf{1}[\text{rank}(t_q) \leq K], \qquad \text{MRR} = \frac{1}{|\mathcal{Q}|}\sum_{q \in \mathcal{Q}} \frac{1}{\text{rank}(t_q)},\] where \(\mathcal{Q}\) is the set of evaluation queries targeting the holdout TFs, and \(\text{rank}(t_q)\) is the rank of the ground-truth tail. Following standard KGC practice, we employ the filtered setting, which excludes other known valid targets for the same query when calculating the rank. This reformulation guarantees that our benchmark directly quantifies a model’s practical capacity to accelerate discovery within realistic experimental budgets.

4 Modeling Gene-Regulation Co‑evolution: The CoDiffGRN Framework↩︎

4.1 Preliminaries↩︎

Diffusion models consist of a forward noising process and a reverse denoising process [27]. The forward process gradually corrupts data \(\mathbf{x}_0\) into noise \(\mathbf{x}_T\) through a Markov chain \(q(\mathbf{x}_{1:T} \mid \mathbf{x}_0) = \prod_{t=1}^T q(\mathbf{x}_t \mid \mathbf{x}_{t-1})\). The reverse Markov process is parameterized by neural networks and progressively denoises the latent states from \(\mathbf{x}_T\) to \(\mathbf{x}_0\). The reverse trajectory is modeled as \[p_\theta(\mathbf{x}_{0:T}) = q(\mathbf{x}_T)\prod_{t=1}^{T} p_\theta(\mathbf{x}_{t-1}\mid \mathbf{x}_t).\]

Diffusion process. We base our formulation on discrete diffusion spaces. Following D3PM [18], we consider one-hot discrete random variables with \(K\) categories, where \(\mathbf{x}_0, \mathbf{x}_1, \dots, \mathbf{x}_t \in \{0,1\}^K\). The forward transition probabilities are encoded by matrices \(Q_t \in \mathbb{R}^{K \times K}\), where \([Q_t]_{ij} = q(x_t = j \mid x_{t-1} = i)\). The transition distribution is defined as: \[q(\mathbf{x}_t \mid \mathbf{x}_{t-1}) = \mathrm{Cat}(\mathbf{x}_t; \mathbf{p} = \mathbf{x}_{t-1} Q_t),\] where \(\mathrm{Cat}(\mathbf{x}; \mathbf{p})\) denotes a categorical distribution over the one-hot vector \(\mathbf{x}\) with probabilities \(\mathbf{p}\). The marginal at timestep \(t\) follows directly from the Markov property: \[q(\mathbf{x}_t \mid \mathbf{x}_0) = \mathrm{Cat}(\mathbf{x}_t; \mathbf{p} = \mathbf{x}_0 \bar{Q}_t), \quad \text{where} \quad \bar{Q}_t = Q_1 Q_2 \cdots Q_t.\]

Reverse process. The reverse process initializes \(\mathbf{x}_T \sim \pi(\mathbf{x})\) and generates \(\mathbf{x}_0\) by reversing the timesteps \(t = T, T-1, \dots, 0\). Following prior work [18], [27], a neural network \(\mathrm{nn}_\theta(\mathbf{x}_t)\) predicts \(\tilde{p}_\theta(\tilde{\mathbf{x}}_0 \mid \mathbf{x}_t)\). The reverse transition is then computed as: \[\label{eq:rev1} p_\theta(\mathbf{x}_{t-1} \mid \mathbf{x}_t) \propto \sum_{\tilde{\mathbf{x}}_0} q(\mathbf{x}_{t-1}, \mathbf{x}_t \mid \tilde{\mathbf{x}}_0) \tilde{p}_\theta(\tilde{\mathbf{x}}_0 \mid \mathbf{x}_t).\tag{1}\] By Bayes’ theorem and the Markov property, the joint posterior admits a tractable form: \[\label{eq:rev2} q(\mathbf{x}_{t-1}, \mathbf{x}_t \mid \mathbf{x}_0) = \frac{q(\mathbf{x}_{t-1} \mid \mathbf{x}_t, \mathbf{x}_0) q(\mathbf{x}_{t-1} \mid \mathbf{x}_0)}{q(\mathbf{x}_t \mid \mathbf{x}_0)} = \mathrm{Cat}\left(\mathbf{x}_{t-1}; \mathbf{p} = \frac{\mathbf{x}_t Q_t^\top \odot \mathbf{x}_0 \bar{Q}_{t-1}}{\mathbf{x}_0 \bar{Q}_t \mathbf{x}_t^\top}\right).\tag{2}\]

Choice of transition matrix and conditioning. Different designs of \(Q_t\) induce distinct inductive biases [18]. A simple yet effective parameterization is \(Q_t = \alpha_t \mathbf{I} + (1 - \alpha_t)\mathbf{1m}^\top\), where \(\mathbf{m}\) denotes the marginal distribution and \(\alpha_t\) controls the noise scale. The cumulative transition admits a closed form \(\bar{Q}_t = \bar \alpha_t \mathbf{I} + (1 - \bar \alpha_t)\mathbf{1m}^\top\), where \(\bar \alpha_t = \prod_{\tau=1}^t \alpha_\tau\). The noise schedule is governed by a cosine function [34]. Furthermore, to enable conditional generation, the reverse process learns a conditional data distribution \(p_\theta(\mathbf{x}_{t-1} \mid \mathbf{x}_t, \mathbf{c})\), where auxiliary condition \(\mathbf{c}\) is injected into the denoiser to guide the reverse trajectory.

4.2 Joint Discrete Diffusion on Gene Regulatory Network↩︎

As illustrated in Figure 2, we propose a joint discrete diffusion framework for GRNs. In contrast to prior works that decouple gene features from edge diffusion [24], [25], our framework unifies discretized gene expression and regulatory interactions into a single co-evolutionary state space.

Cell-cluster discretization. To mitigate feature collapse on unseen genes, we transform continuous scRNA-seq profiles into biologically coherent discrete states. Prior studies show that discretized gene activation patterns capture stable functional structures across diverse biological settings [35], [36]. We cluster cells into \(k\) groups and binarize gene expression within each cluster via thresholding, producing a \(k\)-bit activation pattern per gene. Each binary pattern is mapped to a categorical variable and represented as a one-hot vector in \(\mathbb{R}^{F_D}\), where \(F_D = 2^k\).

Joint modeling. We follow a joint diffusion formulation [37] to unify gene expression and regulatory edges. Edges are represented as \(\mathbf{E} \in \mathbb{R}^{N \times N \times 2}\) (indicating absence/presence), and node features as \(\mathbf{D} \in \mathbb{R}^{N \times F_D}\). The global graph token for each node \(i\) is constructed via concatenation: \(\mathbf{G}[i] = \big[ \mathbf{d}_i \;\Vert\; \mathbf{e}_{i1} \;\Vert\; \cdots \;\Vert\; \mathbf{e}_{iN} \big]\), where \(\mathbf{d}_i \in \mathbb{R}^{F_D}\) denotes the discretized node features, and \(\mathbf{e}_{ij} \in \mathbb{R}^2\) is the one-hot encoding of the edge from node \(i\) to node \(j\).

Design of the transition matrix. We formulate a joint transition matrix \(Q_G \in \mathbb{R}^{F_G \times F_G}\), constructed from four component transition probability matrices: \(Q_D\) (gene \(\to\) gene), \(Q_E\) (regulation \(\to\) regulation), \(Q_{ED}\) (regulation \(\to\) gene), and \(Q_{DE}\) (gene \(\to\) regulation): \[Q_G = \begin{bmatrix} Q_D & \mathbf{1}_N^\top \otimes Q_{DE} \\ \mathbf{1}_N \otimes Q_{ED} & \mathbf{1}_{N \times N} \otimes Q_E \end{bmatrix},\] where \(\otimes\) denotes the Kronecker product, and \(\mathbf{1}_N\), \(\mathbf{1}_{N \times N}\) are all-ones vectors and matrices of appropriate dimensions. Following [38], \(Q_D\) and \(Q_E\) are derived from the marginal distributions \(\mathbf{m}_D\) and \(\mathbf{m}_E\) computed from the training data. For the cross-dependencies, \(\mathbf{m}_{ED}\) and \(\mathbf{m}_{DE}\) capture the conditional co-occurrence probabilities between gene expression types and regulation types. This explicitly ensures that edge transitions are biologically conditioned on the states of their endpoint genes.

Forcing true D technique. During training, the reverse process approximates \(p_\theta(G_{t-1} \mid G_t) \propto \sum_{\tilde{G}_0} q(G_{t-1}, G_t \mid \tilde{G}_0) \tilde{p}_\theta(\tilde{G}_0 \mid G_t)\). To mitigate error accumulation and sharpen the learning signal for unknown edge structures, we apply a teacher-forcing-style technique that replaces the predicted node component \(\tilde{\mathbf{D}}_0\) within the generated \(\tilde{G}_0\) with the ground-truth \(\mathbf{D}_0\).

4.3 Conditional Graph Denoising Architecture↩︎

Conditioning mechanism. We incorporate multiple sources of conditioning information into the denoiser to guide the reverse process. Combining sinusoidal timestep embeddings \(\mathbf{t}_i\), continuous gene expression features \(\mathbf{x}_i\), noisy discrete node states \(\mathbf{d}_{t,i}\), flattened incoming/outgoing edge features \(\mathbf{e}_{t,i:}\), and a boolean TF indicator \(\mathbf{z}_i^{\text{TF}}\), the conditioned input feature for node \(i\) is constructed as: \[\mathbf{h}_i^{(0)} = \big[ \mathbf{t}_i \;\Vert\; \mathbf{x}_i \;\Vert\; \mathbf{d}_{t,i} \;\Vert\; \mathbf{e}_{t,i:} \;\Vert\; \mathbf{z}_i^{\text{TF}} \big].\]

Denoiser backbone. We parameterize the reverse process using a TF-aware graph denoiser built upon GATv2 [39]. Each layer updates node representations dynamically: \[\mathbf{h}_i^{(l+1)} = \sigma\left( \sum_{j \in \mathcal{N}(i)} \alpha_{ij}^{(l)} \mathbf{W}^{(l)} \mathbf{h}_j^{(l)} \right),\] where the attention coefficients \(\alpha_{ij}^{(l)}\) are computed via: \[\alpha_{ij}^{(l)} = \mathrm{softmax}_j \left( \mathbf{a}^\top \, \mathrm{LeakyReLU}\big( \mathbf{W}^{(l)}[\mathbf{h}_i^{(l)} \Vert \mathbf{h}_j^{(l)}] \big) \right).\] The node representations are then projected into TF-specific (\(\mathbf{h}^{\text{tf}}_i\)) and target-specific (\(\mathbf{h}^{\text{tg}}_j\)) embedding spaces. A score head yields a scalar compatibility score \(s_{ij} = f_\mathrm{score}([\mathbf{h}^{\text{tf}}_i \Vert \mathbf{h}^{\text{tg}}_j])\). Binary edge logits are formed as \([-s_{ij},\, s_{ij}]\) and added residually to the noisy edge input \(\mathbf{E}_t\). Node logits are obtained via a linear projection of averaged embeddings with a residual connection to \(\mathbf{D}_t\).

Training objective. The denoiser is trained using a weighted combination of discrete edge reconstruction loss (\(\mathcal{L}_{\text{edge}}\)) and a node-level auxiliary regression objective (\(\mathcal{L}_{\text{node}}\)). Let \(\hat{{e}}_{i,j}\) denote the predicted probability of edge existence and \({e}_{i,j} \in \{0,1\}\) the ground-truth. The total loss is defined as: \[\mathcal{L} = - \sum_{i,j} \left( e_{i,j} \log \hat{e}_{i,j} + (1 - e_{i,j}) \log (1 - \hat{e}_{i,j}) \right) + \lambda \frac{1}{N} \sum_{i=1}^{N} \left\| \hat{\mathbf{X}}_{0,i} - \mathbf{X}_i \right\|_2^2,\] where the auxiliary MSE loss preserves continuous biological signals within the discretized space, and \(\lambda \ll 1\) controls its strength.

4.4 TF-ALL Subgraph Sampling (TASS)↩︎

Diffusion models typically require large training datasets, but GRN inference is inherently data-scarce since each dataset constitutes a single large graph. To address this, we propose TF-ALL Subgraph Sampling (TASS), which constructs a training set by sampling subgraphs from the input graph (Figure 1 (c)). The model is trained on these subgraphs and, at inference time, repeatedly samples subgraphs, predicts edges, and aggregates results to reconstruct the full graph. This expands effective supervision by exposing diverse gene–topology patterns, improving generalization under limited data.

Sampling. As shown in Figure 2, we adopt node-based sampling due to the lack of reliable edge priors, following uniform node sampling [40]. Exploiting the bipartite-like structure of GRNs (TFs \(\to\) targets), each subgraph of size \(k\) is formed by sampling a fraction \(t\) (TF-to-ALL ratio) of nodes from the TF set, and the remainder from all genes, with resampling for TF duplicates. Let \(p=k/n\). Following theoretical bounds for graph recovery [40], we set the number of subgraphs to \(m=\lceil p^{-2}\log n \log(1/\delta)\rceil\). In GRNs, this bipartite sampling strategy inherently yields higher edge coverage due to the sparse, hierarchical nature of regulatory interactions.

Reconstruction. At inference, we apply the same sampling strategy, generate predictions across subgraphs, and map them back to the global edge indices. Outputs are converted to probabilities via softmax, and overlapping predictions are averaged to form the final consensus. Uncovered edges are assigned a default baseline prior.

5 Experiments↩︎

5.1 Experimental Setup↩︎

Evaluation protocols. We primarily evaluate performance using Hits@K and MRR on Inductive Gene-Holdout Split introduced in Section 3. To reflect practical experimental settings, we report Hits@10 and Hits@50, corresponding to low- and high-throughput validation regimes.

Baseline models. We compare against a comprehensive set of baselines, including two traditional methods (GENIE3 [22], GRNBoost2 [20]), two feature-based deep learning models (CNNC [21], DeepSEM [23]), four graph-based models (GNNLink [9], GENELink+ [10], GRNFormer [11], GCLink [12]), and two diffusion-based methods (RegDiffusion [24], DigNet [25]).

Model setups. For fair comparison, all baselines are implemented following their original configurations. Our model adopts the TF-Aware Graph Denoiser, built upon a multi-layer GATv2 architecture with three attention heads per layer and progressively reduced hidden dimensions (128, 64, 64, 32), producing a 16-dimensional output embedding. We use a LeakyReLU with negative slope 0.2 and concatenate multi-head outputs. Score head \(f_\mathrm{score}\) is implemented as a two-layer MLP.

Training setups. All baselines are trained according to their original protocols. Our model is optimized using AdamW with a learning rate of \(1 \times 10^{-4}\), weight decay \(1 \times 10^{-3}\), batch size 64, and trained for 500 epochs. We use a discrete state space with \(k=4\) clusters (yielding 16 classes) and 500 diffusion steps. Subgraphs are sampled with 100 nodes each. All experiments are conducted on a local server with 100 CPU cores and 4 NVIDIA RTX 4090 GPUs (24 GB), and are run over three random seeds (0, 1, and 2), with both mean and standard deviation reported.

3pt

5.2 Results↩︎

We report the performance of representative models on the Specific network under the proposed BEELINE-KGC in Table ¿tbl:tab:benchmark?. Based on these results, we highlight the following key findings:

• The proposed BEELINE-KGC benchmark is substantially more challenging than existing protocols. Traditional and feature-based methods largely collapse under inductive gene holdout, yielding near-zero Hits@10 and MRR. Even several GNN-based models degrade substantially, indicating that BEELINE-KGC reveals inductive generalization failures obscured by transductive benchmarks and global metrics.

• GNN-based models show limited inductive generalization. While GENELink+ remains competitive, most GNN methods suffer from low performance. Their heavy reliance on predefined prior graphs leads to disconnected components for unseen genes, causing error accumulation during message passing and widening the gap against diffusion-based approaches.

• Diffusion-based models offer partial robustness yet remain insufficient. According to Table ¿tbl:tab:benchmark?, diffusion-based methods generally outperform GNN baselines and match GENELink+ in some cases, but still fall behind CoDiffGRN. This suggests that existing diffusion models alleviate reliance on explicit graph structure, yet remain vulnerable to feature collapse on unseen genes. Moreover, their single-graph optimization paradigm appears suboptimal for diffusion.

• CoDiffGRN achieves consistent SOTA performance. CoDiffGRN ranks first across all cell types and metrics in both TFs+500 and TFs+1000 settings, with relative gains up to +184.8% (and +24.5% on average) over the strongest baselines on the most competitive Specific network. These results validate the effectiveness of our co-evolutionary discrete diffusion framework with gene-state-dependent transitions and the TASS strategy for addressing inductive GRN inference and top-\(K\) ranking quality.

5.3 Ablation Study↩︎

To verify the effectiveness of our designs, we conduct extensive ablation experiments on BEELINE-KGC, focusing on the joint diffusion mechanism and the TASS augmentation strategy. The variants are: (i) Edge Only: a variant that only using edge diffusion; (ii) w/o TASS: the model trained on the original graph without TASS; and (iii) RandSS: replacing TASS with uniform random subgraph sampling.

Effect of joint modeling. As shown in Figure 3, removing joint modeling (Edge Only) consistently degrades performance across datasets and metrics. This indicates that decoupling gene expression states from regulatory edges limits the model’s ability to capture conditional regulatory interactions and reduces top-\(K\) ranking quality, especially when generalizing to unseen genes in inductive settings.

Effect of subgraph sampling. Removing TASS leads to a substantial performance drop across all experimental settings. In contrast, uniform random sampling (RandSS) partially recovers the performance and achieves superior results in a few isolated cases, but remains consistently less effective overall. These results highlight the importance of informative subgraph construction and demonstrate that the task-aware sampling strategy of TASS provides more reliable improvements than generic random sampling under inductive and data-scarce GRN inference settings.

6 Conclusion↩︎

In this work, we address the gap between computational GRN inference and real-world biological discovery. We introduce BEELINE-KGC, an application-aligned framework that enforces an inductive, KGC-based evaluation protocol to reflect the top-\(K\) constraints of wet-lab validation. To overcome the exposed gene-level generalization gap, we propose CoDiffGRN, the first joint node-edge discrete diffusion model for GRNs. By capturing the co-evolution of gene expression and network topology, and further enhanced by our TASS strategy, CoDiffGRN achieves robust inductive generalization. Together, our benchmark and model establish a unified paradigm for reliable and budget-aware biological discovery in systems biology.

Limitations and Future Work. Despite its performance, our approach has several limitations that suggest directions for future research. The current evaluation is restricted to established scRNA-seq systems, which may not fully represent the complexity of rare tissues or cross-species regulatory dynamics. We plan to extend our benchmark to include more diverse biological contexts and validate the model’s transferability across different species. Moreover, the joint diffusion process and TASS strategy introduce non-trivial computational overhead during training and inference. Future work will focus on developing optimized sampling kernels and parallelized subgraph processing to enhance scalability.

7 Results on Previous Transductive Setting↩︎

3pt

3pt

8 Cell-cluster Discretization↩︎

This section provides additional details and supporting evidence for the proposed Cell-cluster-based Discretization strategy. We first present the algorithmic procedure, followed by visualizations that demonstrate the biological coherence of the induced discrete gene state space.

8.1 Algorithm↩︎

8.2 Design Motivation and Biological Evidence↩︎

We provide the design motivation and biological evidence supporting the proposed Cell-cluster-based Discretization. The goal of this design is to construct a biologically coherent discrete state space that enables stable inductive diffusion over unseen genes.

Motivation. A natural alternative is to directly cluster gene expression profiles using standard unsupervised methods (e.g., k-means). However, such approaches may produce clusters dominated by statistical artifacts such as sparsity patterns or zero inflation in scRNA-seq data, rather than biologically meaningful structure. In contrast, scRNA-seq measurements are inherently organized by cellular contexts (e.g., conditions or states), suggesting that discretization should be grounded in cell-level structure to preserve functional coherence.

Design Principle. Our method explicitly leverages this structure by first clustering cells into \(k\) biologically meaningful groups, and then constructing gene-level discrete representations based on expression consistency across these cell states. This ensures that discretization reflects stable functional activity rather than superficial similarity in gene expression vectors.

Figure Organization and Evidence. We use UMAP [41] to visualize the resulting cell clusters obtained from scRNA-seq data in Figure 5. The clusters correspond to distinct cellular states, indicating that the clustering captures meaningful biological organization rather than arbitrary partitions.

Figure 6 presents representative genes across different cell clusters. We observe that genes exhibit consistent activation patterns within clusters, suggesting that discretization captures stable gene-level functional behavior conditioned on cellular context. This is further proved by Figure 7, which visualizes the induced gene representation manifold in the discrete state space using UMAP [41]. The resulting structure exhibits smooth organization and local coherence, suggesting that the discretization preserves biologically meaningful continuity while providing a structured state space suitable for diffusion modeling.

Together, these results support the role of Cell-cluster-based Discretization as a biologically grounded representation scheme that enables stable inductive generalization under sparse and heterogeneous scRNA-seq data.

9 TF-ALL Subgraph Sampling (TASS)↩︎

In this section, we present the complete pipeline of TF-ALL Subgraph Sampling (TASS), which consists of a stochastic subgraph sampling strategy followed by a consensus-based reconstruction procedure for global edge inference.

9.1 Sampling↩︎

9.2 Reconstruction↩︎

10 Dataset Information of BEELINE-KGC↩︎

The original BEELINE benchmark is constructed from five experimental single-cell RNA-seq datasets covering seven cell types (hESC, hHEP, mDC, mESC, mHSC-E, mHSC-GM, mHSC-L) across human and mouse systems [14]. These datasets span diverse biological contexts, including differentiation and developmental processes, and provide measurements that enable pseudotime inference.

Each dataset is paired with up to four types of ground-truth networks:: cell-type-specific (Specific) network derived from ChIP-seq experiments and curated databases such as ENCODE [42], ChIP-Atlas [43], and ESCAPE [44]; non-specific (Non-Specific) network aggregated from large-scale resources such as DoRothEA [45], RegNetwork [46], and TRRUST [47]; functional interaction networks from STRING [48], capturing broader biological associations that may include indirect regulatory relationships; loss-/gain-of-function (LOF/GOF) evidence from ESCAPE [44] (available only for mESC). For each dataset, two gene subsets are constructed by selecting the top 500 or 1,000 most variably expressed genes with transcription factors (Abbr.: TFs+500 and TFs+1000, respectively), culminating in total 44 distinct evaluation tasks.

Among these settings, future work may place particular emphasis on the Specific networks, as their cell-type-matched regulatory evidence provides the closest alignment between the expression measurements and the evaluation targets, making them a more representative testbed for context-dependent GRN discovery. The other network types remain valuable as complementary settings for assessing robustness to broader, heterogeneous, or functional interaction evidence.

References↩︎

---

1. Corresponding author↩︎