Paper deep dive
$D^{2}R^{2}$: Discrete Diffusion with Regulation Reinforcement for Single-Cell Perturbation Prediction
Ninghan Fan, Qi Liu, Xunuo Zhu, Yukai Sun, Luyuan Chen, Xuheng Zhou, Yuetian Du, Ming Kong, Xiaojun Zhu, Jie Liu, Zhan Zhou, Qiang Zhu
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 8/18/2026, 6:04:57 AM
Summary
The paper introduces D2R2, a novel framework for single-cell perturbation prediction that reformulates the task as regulation-guided gene-wise progressive generation. Unlike existing methods that predict entire expression profiles simultaneously, D2R2 uses a Masked Discrete Diffusion Model (MDDM) to reconstruct expression tokens step-by-step. A Regulatory Policy Module (RPM) determines the order in which genes are generated, initialized from a gene regulatory network (GRN) and refined using Group-Relative Policy Optimization (GRPO). This approach leverages biological priors to prioritize regulatory genes early, improving prediction accuracy and interpretability on benchmarks like Norman19 and VCC-H1.
Entities (8)
Relation Signals (7)
D2R2 → evaluatedon → Norman19
confidence 96% · Across Norman19 and VCC-H1, D2R2 achieves the best performance on all five metrics on Norman19
D2R2 → evaluatedon → VCC-H1
confidence 96% · Across Norman19 and VCC-H1, D2R2 achieves the best performance on all five metrics on Norman19 and remains competitive on H1.
D2R2 → uses → Regulatory Policy Module
confidence 95% · D2R2 couples a Masked Discrete Diffusion Model (MDDM) for progressive expression generation with a Regulatory Policy Module (RPM)
D2R2 → uses → Masked Discrete Diffusion Model
confidence 95% · D2R2 couples a Masked Discrete Diffusion Model (MDDM) for progressive expression generation with a Regulatory Policy Module
Regulatory Policy Module → optimizesusing → Group Relative Policy Optimization
confidence 92% · Then, group-relative policy optimization refines only the ordering policy using final perturbation-effect agreement as reward.
Regulatory Policy Module → initializesfrom → Gene Regulatory Network
confidence 90% · A Regulatory Policy Module initializes the generation policy from a gene regulatory network inferred from control cells
Gene Regulatory Network → inferredusing → DeepSEM
confidence 85% · we infer a cell-type-specific GRN Gz=(S,Az) from unperturbed control cells using DeepSEM
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Predicting single-cell transcriptomic responses to genetic perturbations is central to functional genomics and virtual-cell modeling. Existing approaches, however, typically predict an entire expression profile as a whole, leaving the order in which individual gene responses are generated unmodeled. To address this problem, we introduce \textbf{$D^{2}R^{2}$} (\textbf{D}iscrete \textbf{D}iffusion with \textbf{R}egulation \textbf{R}einforcement), which reformulates perturbation prediction as regulation-guided gene-wise progressive generation. A Masked Discrete Diffusion Model represents expression as ordinal tokens and reconstructs a fully masked profile step by step, allowing generated gene responses to condition those that remain masked. A Regulatory Policy Module initializes the generation policy from a gene regulatory network inferred from control cells and adapts it to the perturbation and current partially generated state. Then, group-relative policy optimization refines only the ordering policy using final perturbation-effect agreement as reward. Across Norman19 and VCC-H1, $D^{2}R^{2}$ achieves the best performance on all five metrics on Norman19 and remains competitive on H1. Controlled ablations holding the generator and generation budget fixed show that biological-prior ordering improves over random ordering and is more reliable than uncertainty-based heuristics, whereas reversing the biological-prior ordering degrades every metric. Biological analyses further show that the refined policy prioritizes regulatory genes early while promoting perturbation-specific transcription factors and responsive genes. These results establish gene generation order as an effective, controllable, and biologically interpretable dimension of single-cell perturbation prediction.
Tags
Links
- Source: https://arxiv.org/abs/2608.15288v1
- Canonical: https://arxiv.org/abs/2608.15288v1
Trouble viewing inline? Open PDF directly →
Full Text
63,703 characters extracted from source content.
Expand or collapse full text
D2R2D^2R^2: Discrete Diffusion with Regulation Reinforcement for Single-Cell Perturbation Prediction Ninghan Fan Qi Liu Xunuo Zhu Yukai Sun Luyuan Chen Xuheng Zhou Yuetian Du Ming Kong Xiaojun Zhu Jie Liu Zhan Zhou Qiang Zhu Abstract Predicting single-cell transcriptomic responses to genetic perturbations is central to functional genomics and virtual-cell modeling. Existing approaches, however, typically predict an entire expression profile as a whole, leaving the order in which individual gene responses are generated unmodeled. To address this problem, we introduce D2R2D^2R^2 (Discrete Diffusion with Regulation Reinforcement), which reformulates perturbation prediction as regulation-guided gene-wise progressive generation. A Masked Discrete Diffusion Model represents expression as ordinal tokens and reconstructs a fully masked profile step by step, allowing generated gene responses to condition those that remain masked. A Regulatory Policy Module initializes the generation policy from a gene regulatory network inferred from control cells and adapts it to the perturbation and current partially generated state. Then, group-relative policy optimization refines only the ordering policy using final perturbation-effect agreement as reward. Across Norman19 and VCC-H1, D2R2D^2R^2 achieves the best performance on all five metrics on Norman19 and remains competitive on H1. Controlled ablations holding the generator and generation budget fixed show that biological-prior ordering improves over random ordering and is more reliable than uncertainty-based heuristics, whereas reversing the biological-prior ordering degrades every metric. Biological analyses further show that the refined policy prioritizes regulatory genes early while promoting perturbation-specific transcription factors and responsive genes. These results establish gene generation order as an effective, controllable, and biologically interpretable dimension of single-cell perturbation prediction. Introduction Figure 1: Paradigm comparison. Each row summarizes paradigm properties and traces gene-expression generation, with web nodes and aligned vector entries denoting genes and their expression values. We compare traditional single-pass prediction, DDPM-based global denoising, and two forms of D2R2D^2R^2: random unmasking (MDDM) and regulation-guided unmasking. Accurate prediction of cellular responses to perturbations is a central challenge in functional genomics and drug discovery. Despite rapid advances in high-throughput perturbational measurements (Qian, Dong, and Guo 2025), experimentally producing perturbation data remains costly and laborious. These challenges have motivated the development of virtual cell models that simulate cellular responses in silico, offering the potential to advance biological understanding and accelerate personalized medicine (Bunne et al. 2024; Roohani et al. 2025). However, accurately modeling perturbation-induced transcriptomic responses remains difficult due to complex gene regulatory mechanisms, biological heterogeneity, and dynamic cellular responses. Methods for perturbation prediction have evolved from single-pass predictors built on transformers (Adduri et al. 2025; Cui et al. 2024b) or variational autoencoders (Bereket and Karaletsos 2023) to diffusion-based generative models (He et al. 2026; Klein et al. 2025; Yu et al. 2026), as illustrated in Fig. 1a and b. Despite their architectural diversity, these approaches predominantly follow a expression profile-level generation paradigm: single-pass models generate the full profile at once, whereas diffusion models iteratively denoise the entire profile. However, in both cases, the gene-wise generation order remains unmodeled. This omission matters because gene responses are not independent, predicting gene responses in a structured order would allow earlier predictions to provide informative context for later ones (Freimer et al. 2022; Song et al. 2025). This motivates a new perturbation-prediction paradigm that explicitly models gene-wise resolution order, allowing individual gene responses to be generated progressively. Decomposing profile generation into gene-wise generation introduces a new challenge: in what order should gene responses be generated? Simple strategies such as random selection, predictive confidence (Chang et al. 2022), or entropy (Settles 2009) use either no structural information or only the generator’s instantaneous uncertainty. A biologically grounded ordering could instead exploit gene regulatory structure to identify genes likely to influence the responses of others. Yet any fixed regulatory prior derived from unperturbed control cells would capture cell-type-level structure rather than perturbation-specific dependencies. The same ordering could not determine which dependencies are most useful for a particular perturbation in the same cell line. Addressing these limitations requires an ordering mechanism that is biologically grounded at initialization yet able to adapt to both the perturbation and the evolving generation state. To overcome the aforementioned challenges, we reformulate single-cell perturbation prediction as an gene-wise progressive generation. Based on this formulation, we propose D2R2D^2R^2, a Discrete Diffusion model with Regulation Reinforcement, which couples a Masked Discrete Diffusion Model (MDDM) for progressive expression generation with a Regulatory Policy Module (RPM) for biologically grounded, adaptive gene ordering, as shown in Fig. 1c and d. Specifically, MDDM represents expression values as ordered discrete tokens and begins generation from a fully masked profile. At each step, RPM selects a subset of unresolved genes, and MDDM predicts their expression tokens. Once generated, these tokens remain fixed and provide explicit context for subsequent steps. RPM determines which genes to generate next by transforming a control-derived regulatory prior into a state- and perturbation-specific generation policy. Training proceeds in two stages. We first randomly mask expression tokens at sampled diffusion times and train MDDM with an x0x_0-prediction objective to reconstruct the original tokens from partially masked profiles. We then freeze MDDM, initialize RPM to reproduce an ordering prior derived from a gene regulatory network inferred from unperturbed control cells, and refine only the ordering policy with group-relative policy optimization (GRPO), using final perturbation-response agreement as the reward. Across Norman19 and VCC-H1, experiments validate the effectiveness of D2R2D^2R^2. Controlled ablations with a fixed MDDM show that generation order is consequential: the biological prior outperforms random and uncertainty-based strategies, whereas reversing it consistently degrades prediction. RPM further adapts this prior to the perturbation and evolving generation state, while biological analyses show that the refined orders preserve early regulatory-gene prioritization and promote perturbation-specific transcription factors (TFs) and responsive genes. Together, these results establish generation order as an effective and biologically interpretable component of perturbation prediction. Our contributions are threefold: • We formulate single-cell perturbation prediction as gene-wise progressive generation with MDDM, making the generation order explicit and allowing already-resolved gene responses to condition those that remain masked. • We introduce RPM, which initializes an ordering policy from a control-derived biological prior and refines it with GRPO while keeping the generator fixed, producing state- and perturbation-specific generation orders. • We conduct a controlled study of random, uncertainty-based, biological, reversed, and refined orders, establishing generation order and its directionality as substantive design choices. Related Work Single-cell perturbation prediction and generative modeling. Single-cell perturbation prediction has been studied under diverse predictive and generative formulations. Latent-variable methods such as scGen model perturbation effects as shifts in a shared cellular space (Lotfollahi, Wolf, and Theis 2019), while disentanglement-based predictors separate cellular state, perturbation, and covariate factors to improve generalization (Piran et al. 2024). Distributional and flow-based approaches instead learn mappings between control and perturbed cell populations (Bunne et al. 2023; Klein et al. 2025; Yu et al. 2026). More recently, foundation models use large-scale single-cell pretraining to obtain transferable representations for perturbation prediction (Cui et al. 2024a; Theodoris et al. 2023; Adduri et al. 2025). Although these approaches differ in their representations and predictive targets, they produce the response through a forward prediction, latent transformation, or distributional map. Diffusion-based models introduce iterative generation by denoising cellular states or learning bridges between control and perturbed populations (Tang et al. 2023; He et al. 2026). Their trajectories, however, are typically governed by noise schedules over the whole transcriptomic state rather than an explicit gene-wise resolution order. D2R2D^2R^2 treats perturbation prediction as masked discrete generation over expression tokens, making this order a controllable part of the model. Biological priors and regulatory dependencies. Cellular responses to perturbations are constrained by biological structure, including gene regulatory networks (GRN), transcription factor-target relationships, and pathway interactions (Dixit et al. 2016; Norman et al. 2019; Freimer et al. 2022; Shu et al. 2021; Yuan et al. 2021). Prior work incorporates such information through graph neural networks, message passing, structure-aware embeddings, or adjacency-based regularization (Adduri et al. 2025; Theodoris et al. 2023); graph-based perturbation models similarly use gene–gene relationships to improve response prediction (Roohani, Huang, and Leskovec 2024). These approaches demonstrate the value of regulatory information, but primarily use it to shape gene representations or information exchange within the predictor. Regulatory structure therefore affects how genes interact during prediction, without defining a measurable sequence in which their responses are resolved. D2R2D^2R^2 applies the structure at this different point: it uses a control-derived GRN to initialize the order in which gene responses are committed, after which the ordering policy is optimized for perturbation prediction. Reinforcement learning for adaptive generation policies. Reinforcement learning is useful in computational biology when local decisions are difficult to supervise but complete designs can be evaluated through a delayed objective. Applications include optimizing molecular properties (Olivecrona et al. 2017), constructing molecular graphs under non-differentiable objectives (You et al. 2018), RNA inverse folding (Runge et al. 2018), and protein design (Lutz et al. 2023). A parallel line applies RL to sequential generation itself, using sequence-level rewards or actor–critic objectives to judge early actions by their downstream consequences (Ranzato et al. 2015; Bahdanau et al. 2016). D2R2D^2R^2 has the same delayed-reward structure, but the policy acts on generation order rather than expression values: its actions select genes to unmask, and their value depends on how well the resulting predictions serve as context for the genes generated later. This allows RPM to optimize the ordering policy while the expression generator remains fixed. Method Figure 2: Overview of the D2R2D^2R^2 framework. The left panel shows the iterative generation loop: after the MDDM predicts the selected expression tokens, RPM observes the updated partial state and recomputes which genes should be resolved next. The upper-right panel shows the biological prior used to initialize the policy, while the lower-right panel illustrates GRPO that optimizes the order to each perturbation. D2R2D^2R^2 first learns to generate perturbed expression tokens with a masked discrete diffusion model, and then learns which unresolved genes should be predicted at each generation step. This separates expression modeling from generation-order learning and allows us to vary the order while holding the generator fixed. Preliminaries Problem formulation. We consider perturbation response prediction over a fixed gene set S=g1,…,gLS=\g_1,…,g_L\ of L genes, indexed by ℐ=[L]=1,…,LI=[L]=\1,…,L\. For each sample, let c∈ℝLc ^L denote the control expression profile, y∈ℝLy ^L denote the ground-truth perturbed expression profile, p denote the perturbation condition, and z denote cellular covariates. The goal is to learn a conditional generator that predicts y from the conditioning tuple u=(c,p,z)u=(c,p,z). Discrete expression tokens. D2R2D^2R^2 represents the perturbed response in a discrete expression-token space, mapping each continuous expression value yiy_i to an ordinal token x0,i=Q(yi)∈0,…,B−1x_0,i=Q(y_i)∈\0,…,B-1\, where Q(⋅)Q(·) is an equal-frequency binning operator with B bins. We refer to x0=(x0,1,…,x0,L)x_0=(x_0,1,…,x_0,L) as the clean discrete perturbed state, i.e., the tokenized form of y, where position i is tied to the fixed gene identity gig_i. A special mask token m denotes an unpredicted expression token. The mask token makes the unresolved/resolved status of each gene explicit, thereby allowing the gene generation order to be directly specified and optimized. Masked Discrete Diffusion Model D2R2D^2R^2 uses a masked discrete diffusion model as the token-level generator (Austin et al. 2021; Nie et al. 2025a; Nie et al. 2025b), naturally allowing already generated genes to serve as explicit context for masked genes. It learns the conditional distribution pθ(x0∣u)p_θ(x_0 u), where θ denotes the generator parameters. Forward masking corruption. During training, we construct a partially masked input xtx_t from the ground-truth token state x0x_0. Following masked diffusion models (Austin et al. 2021; Nie et al. 2025a), each expression token is independently retained with probability αt=1−t _t=1-t and replaced with the mask token m otherwise: qt(xt∣x0)=∏i=1L[αtδx0,i(xt,i)+(1−αt)δm(xt,i)].q_t(x_t x_0)= _i=1^L [ _t _x_0,i(x_t,i)+(1- _t) _m(x_t,i) ]. (1) Here, δa(x) _a(x) denotes one-hot distribution at a, taking value 11 when x=ax=a and 00 otherwise. This forward process enables continuous-time data corruption, with each gene xt,ix_t,i at position i equaling the original token x0,ix_0,i with probability αt _t and the mask token m with probability 1−αt1- _t. Conditional denoising parameterization. The generator fθf_θ is a non-causal Transformer (Vaswani et al. 2017) that predicts the original token at each masked gene position. We encode the control profile, perturbation condition, and cellular covariates as Hc=Encc(c)H_c=Enc_c(c), Hp=Embp(p)H_p=Emb_p(p), and Hz=Encz(z)H_z=Enc_z(z), and concatenate them into Hcond=[Hc;Hp;Hz]H_cond=[H_c;H_p;H_z]. For each fixed gene position i∈ℐi , the model input is Ht,i=Embx(xt,i)+Embg(gi)H_t,i=Emb_x(x_t,i)+Emb_g(g_i). For each bin b∈0,…,B−1b∈\0,…,B-1\, the predictor outputs pθ(x0,i=b∣xt,u)=Softmaxb(fθ(Ht,Hcond)i).p_θ(x_0,i=b x_t,u)=Softmax_b\! (f_θ(H_t;H_cond)_i ). (2) Training objective. We train fθf_θ with an x0x_0-prediction objective on masked positions. Given x0x_0, we sample t∼Unif(0,1)t (0,1) and xt∼qt(xt∣x0)x_t q_t(x_t x_0), and let ℳt=i:xt,i=mM_t=\i:x_t,i=m\ denote the masked gene positions. Writing pmask(t)=1−αtp_mask(t)=1- _t, the masked diffusion loss is ℒMDDM(θ)=x0,t,xt[ℓt], _MDDM(θ)=E_x_0,t,x_t\! [ _t ], (3) ℓt=1pmask(t)∑i∈ℳt−logpθ(x0,i∣xt,u). _t= 1p_mask(t) _i _t- p_θ(x_0,i x_t,u). In practice, it is estimated by masked cross-entropy over ℳtM_t (Nie et al. 2025a; Nie et al. 2025b). Iterative denoising. At inference time, D2R2D^2R^2 starts from the fully masked state x1=(m,…,m)x_1=(m,…,m): at each step the predictor outputs token distributions for the currently masked positions, selected positions are updated with sampled tokens, and once a token is predicted it is kept fixed. After generation, the discrete token sequence is converted into a continuous expression profile through a deterministic mapping that assigns each token the representative value associated with its bin. The MDDM therefore determines the expression-token distribution for every unresolved gene, but which genes should be committed next iteration is still an open question, leaving space for exploitation. Gene Generation Order Generation-order strategies. We formulate this open choice as a generation-order strategy that determines which masked genes are committed at each iteration. An ordering strategy assigns a priority to every position in the current masked set ℳkM_k and selects a subset ak⊆ℳka_k _k to resolve. This formulation accommodates simple orderings such as Random, as well as model-uncertainty strategies in Table 2: Confidence high/low ranks genes by the largest predicted token probability, while Entropy high/low ranks them by predictive entropy; “high” and “low” select the largest and smallest scores, respectively (Chang et al. 2022; Settles 2009). These strategies provide useful controls, but they use either no structural information or only the generator’s instantaneous uncertainty; neither represents how a perturbation response may propagate across genes. Biological-prior ordering. To introduce regulatory structure into the generation process, we infer a cell-type-specific GRN z=(S,Az)G_z=(S,A_z) from unperturbed control cells using DeepSEM (Shu et al. 2021). The weighted adjacency matrix AzA_z represents directed regulatory interactions among genes. We apply PageRank to this weighted graph to obtain a regulatory-influence score for each gene, which accounts for both its direct connections and its recursive influence through the GRN. Sorting genes by decreasing PageRank score yields a fixed biological-prior order, with influential upstream regulators placed earlier. At generation step k, this prior is restricted to the unresolved genes in ℳkM_k to define the biological supervision qkbioq_k^bio used to initialize RPM. The resulting policy encourages regulatory genes to be generated first, so that their predicted responses provide context for genes resolved later. Perturbation-specific regulatory policy. We introduce the Regulatory Policy Module (RPM) to transform the static cell-type prior into a state- and perturbation-specific generation policy. At step k, the partially generated state is xkx_k, the unresolved positions are ℳk=i:xk,i=mM_k=\i:x_k,i=m\, and the policy state is sk=(xk,p,z,ℳk)s_k=(x_k,p,z,M_k). RPM encodes the current state together with the perturbation and cellular context and produces an unmasking probability for every unresolved gene: πψ(i∣sk)=exp(hψ,i(sk)/τ)∑j∈ℳkexp(hψ,j(sk)/τ),i∈ℳk, _ψ(i s_k)= (h_ψ,i(s_k)/τ) _j _k (h_ψ,j(s_k)/τ), i _k, (4) where hψ,i(sk)h_ψ,i(s_k) is the RPM priority logit and τ is a temperature. The subset size follows a predefined schedule with ∑k|ak|=L _k|a_k|=L. Crucially, RPM recomputes the distribution after every update: once aka_k is filled by the MDDM, the resulting state xk+1x_k+1 becomes the context for deciding which genes to resolve next. The architecture used to parameterize hψh_ψ is described in the Appendix. We summarize the whole process in Algorithm 1. 1 Input: u=(c,p,z)u=(c,p,z), MDDM pθp_θ, RPM πψ _ψ, and bkk=1K\b_k\_k=1^K satisfying ∑k=1Kbk=L _k=1^Kb_k=L. 2 3 Initialize: x1=(m,…,m)x_1=(m,…,m) and ℳ1=ℐM_1=I. 4 5 for k=1,…,Kk=1,…,K do 6 7 8 9 1. Generate gene tokens 10 11 x^k,i∼pθ(x0,i∣xk,u),i∈ℳk x_k,i p_θ(x_0,i x_k,u), i _k. 12 13 14 2. Select gene tokens 15 16 ak←TopKbkπψ(i∣xk,p,z,ℳk):i∈ℳka_k _b_k \ _ψ(i x_k,p,z,M_k):i _k \. 17 18 19 3. Update 20 21 xk+1←xk,xk+1,ak←x^k,akx_k+1← x_k, x_k+1,a_k← x_k,a_k. 22 23 ℳk+1←ℳk∖akM_k+1 _k a_k. 24 25 26 end for 27 28 Return: bin-representative decoding of xK+1x_K+1. 29 30 Algorithm 1 Step-adaptive generation with RPM Biological initialization and GRPO refinement. RPM is first initialized to reproduce the biological-prior policy. For policy states encountered during generation, we minimize ℒbio(ψ)=−sk[∑i∈ℳkqk,ibiologπψ(i∣sk)],L_bio(ψ)=-E_s_k [ _i _kq_k,i^bio _ψ(i s_k) ], (5) which supplies a stable starting point rather than learning from scratch. We then freeze the MDDM and biological prior and refine RPM with GRPO (Shao et al. 2024). Next, we construct G branched rollouts from shared generation state s. In rollout r, RPM samples an ordered subset a(r)=(i1,…,ib)a^(r)=(i_1,…,i_b) without replacement, with likelihood Pψ(a(r)∣s)=∏j=1bπψ(ij∣s)1−∑ℓ<jπψ(iℓ∣s).P_ψ(a^(r) s)= _j=1^b _ψ(i_j s)1- _ <j _ψ(i_ s). (6) The sampled genes are resolved by MDDM, and the remaining trajectory is completed by the frozen reference policy πref _ref. Each rollout receives R(r)=corr(y^(r)−c,y−c)R^(r)=corr( y^(r)-c,y-c), which is normalized within the group to obtain A(r)A^(r). Let ρ(r)=Pψ(a(r)∣s)/Pψold(a(r)∣s)ρ^(r)=P_ψ(a^(r) s)/P_ _old(a^(r) s) and ρ¯(r)=clip(ρ(r),1−ϵ,1+ϵ) ρ^(r)=clip(ρ^(r),1-ε,1+ε). RPM is optimized using ℒGRPO= _GRPO= −r[min(ρ(r)A(r),ρ¯(r)A(r))] -E_r [ (ρ^(r)A^(r), ρ^(r)A^(r) ) ] (7) +βDKL(πψ∥πref). +β D_KL ( _ψ\,\|\, _ref ). Training samples subsets without replacement, whereas inference recomputes RPM priorities at each step and applies deterministic Top-K selection. Experiments Experimental Setup Datasets. We evaluate D2R2D^2R^2 on Norman19 and VCC-H1(H1). Norman19 is a Perturb-seq dataset of gene overexpression perturbations in K562 cells, containing 287 perturbation conditions, including 131 combinatorial perturbations (Norman et al. 2019). VCC-H1 is the H1 human embryonic stem cell benchmark from the Virtual Cell Challenge, containing approximately 400,000 single cells across 300 CRISPRi perturbation targets (Roohani et al. 2025). Norman19 uses a held-out-combination split to test generalization to unseen perturbation combinations, whereas H1 uses a held-out-cell split to evaluate large-scale response reconstruction under the observed perturbation distribution. The two benchmarks therefore cover complementary evaluation regimes. Baselines. Generative models form the primary comparison group: SAMS-VAE (Bereket and Karaletsos 2023), Cell Flow (Klein et al. 2025), Squidiff (He et al. 2026), and scDFM (Yu et al. 2026), covering VAE-based, flow-based, diffusion-based, and distribution-level approaches. We also compare with the perturbation-response predictor Biolord (Piran et al. 2024) and the virtual-cell model STATE (Adduri et al. 2025). All methods use matched splits, preprocessing, and evaluation; implementation details are provided in the Appendix. Metrics. Following PerturBench (Wu et al. 2024), we evaluate differential expression using Pearson Δ , Cos. LogFC, and Cos. LogFC Rank. These metrics measure expression changes relative to control cells, testing whether the model recovers the transcriptional programs activated or suppressed by a perturbation; Cos. LogFC Rank additionally measures whether responses to different perturbations remain distinguishable. For absolute expression, Cos. PCA measures centroid agreement between predicted and observed cells, testing whether the overall post-perturbation cellular state is reconstructed. Sym. KL further evaluates agreement between cell-population distributions, capturing response heterogeneity beyond the centroid. Metric details are provided in the Appendix. Main Results Table 1: Main results on Norman19 and H1. All methods are evaluated over five independent training runs using seeds 0,1,2,3,42\0,1,2,3,42\, and results are reported as mean ± standard deviation. For D2R2D^2R^2, MDDM checkpoint is fixed across runs, while RPM is independently trained for each seed. Metrics are grouped by differential and absolute expression. Best and second-best results are bolded and underlined, respectively. Differential Expression Absolute Expression Dataset Method Pearson Δ ↑ Cos. LogFC ↑ Cos. LogFC Rank ↓ Cos. PCA ↑ Sym. KL ↓ Norman19 SAMS-VAE (NeurIPS 2023) 0.553± 0.079 0.560± 0.043 0.028± 0.010 0.251± 0.042 50.888± 4.215 Biolord (Nat. Biotechnol. 2024) 0.404± 0.069 0.388± 0.051 0.046± 0.010 0.035± 0.016 0.526± 0.040 STATE (NeurIPS 2025) 0.467± 0.020 0.379± 0.020 0.525± 0.017 -0.009± 0.039 0.273± 0.014 Cell Flow (bioRxiv 2025) 0.489± 0.077 0.373± 0.066 0.329± 0.013 0.084± 0.066 0.207± 0.024 Squidiff (Nat. Methods 2025) 0.291± 0.001 0.226± 0.002 0.427± 0.007 0.110± 0.000 2.861± 0.019 scDFM (ICLR 2026) 0.526± 0.010 0.313± 0.015 0.379± 0.008 0.264± 0.032 0.390± 0.034 D2R2D^2R^2 (Ours) 0.706± 0.004 0.724± 0.006 0.024± 0.006 0.703± 0.009 0.148± 0.007 H1 SAMS-VAE (NeurIPS 2023) 0.262± 0.017 0.231± 0.016 0.192± 0.028 0.319± 0.042 40.638± 3.679 Biolord (Nat. Biotechnol. 2024) 0.201± 0.008 0.159± 0.008 0.063± 0.007 0.151± 0.010 0.240± 0.018 STATE (NeurIPS 2025) 0.253± 0.012 0.178± 0.006 0.095± 0.006 0.378± 0.029 0.272± 0.016 Cell Flow (bioRxiv 2025) 0.042± 0.016 0.028± 0.009 0.493± 0.013 0.014± 0.009 0.089± 0.006 Squidiff (Nat. Methods 2025) 0.046± 0.000 0.049± 0.012 0.368± 0.032 -0.008± 0.005 8.435± 0.372 scDFM (ICLR 2026) 0.069± 0.004 0.038± 0.009 0.421± 0.029 0.029± 0.010 0.390± 0.025 D2R2D^2R^2 (Ours) 0.271± 0.002 0.541± 0.004 0.001± 0.000 0.358± 0.001 0.229± 0.002 Table 1 shows that D2R2D^2R^2 ranks first across all metrics on Norman19 and leads in most metrics on H1. In particular, it leads all three differential-expression metrics on both benchmarks, demonstrating that it accurately captures perturbation-induced transcriptional changes while preserving the distinctions among responses to different perturbations. For absolute expression, D2R2D^2R^2 achieves the best centroid agreement on Norman19 and remains competitive on H1, indicating effective reconstruction of the overall post-perturbation cellular state. Its strong Sym. KL results further show that this agreement extends beyond the centroid to the cell-population distribution. Beyond the aggregate results, Fig. 3 compares the per-perturbation distributions of Pearson Δ and Cos. LogFC on Norman19. D2R2D^2R^2 shifts both distributions toward higher response similarity, indicating that its gains are consistent across perturbations. Figure 3: Violin plots showing the distribution of evaluation metrics across perturbations in Norman19. Mann-Whitney U test, * P < 0.05, *** P < 0.001. Table 2: Generation-order ablation on Norman19. All methods use the same frozen MDDM. Random averages multiple random orders. Cells show score and relative performance gain over Random, with signs aligned so that positive always indicates improvement; for Cos. LogFC Rank and Sym. KL, this corresponds to relative error reduction. Orange intensity indicates a positive gain. Percentages are computed before rounding. Order Pearson Δ ↑ Cos. LogFC ↑ Cos. LogFC Rank ↓ Cos. PCA ↑ Sym. KL ↓ Random 0.623 (ref.) 0.720 (ref.) 0.032 (ref.) 0.597 (ref.) 0.183 (ref.) Model-Uncertainty Orders Confidence high 0.590 (−5.3-5.3%) 0.537 (−25.4-25.4%) 0.030 (+5.2+5.2%) 0.353 (−40.9-40.9%) 0.410 (−124.4-124.4%) Confidence low 0.672 (+7.9+7.9%) 0.522 (−27.5-27.5%) 0.068 (−114.1-114.1%) 0.640 (+7.2+7.2%) 0.203 (−11.0-11.0%) Entropy high 0.624 (+0.2+0.2%) 0.452 (−37.2-37.2%) 0.073 (−130.9-130.9%) 0.564 (−5.6-5.6%) 0.206 (−12.5-12.5%) Entropy low 0.591 (−5.1-5.1%) 0.538 (−25.3-25.3%) 0.030 (+5.2+5.2%) 0.355 (−40.5-40.5%) 0.399 (−118.5-118.5%) Biological-Prior Orders Biological prior 0.706 (+13.4+13.4%) 0.703 (−2.3-2.3%) 0.027 (+14.4+14.4%) 0.699 (+17.0+17.0%) 0.119 (+35.0+35.0%) Reversed prior 0.573 (−8.1-8.1%) 0.533 (−26.0-26.0%) 0.064 (−101.8-101.8%) 0.509 (−14.7-14.7%) 0.197 (−7.5-7.5%) Adaptive RPM Policies (GRPO) RPM (random init.) 0.677 (+8.6+8.6%) 0.734 (+2.0+2.0%) 0.026 (+19.0+19.0%) 0.675 (+13.1+13.1%) 0.186 (−1.6-1.6%) RPM (bio init.) 0.711 (+14.2+14.2%) 0.723 (+0.5+0.5%) 0.020 (+35.8+35.8%) 0.693 (+16.1+16.1%) 0.138 (+24.5+24.5%) Generation-Order Ablation To isolate the effect of generation order, we compare different orders. Table 2 reports a controlled ablation in which every configuration uses the same frozen MDDM and generation budget; the only difference is which masked genes are selected at each step. Model uncertainty does not identify a useful generation order. Overall, model-uncertainty-based strategies provide little guidance for organizing gene generation. Although some uncertainty orders improve individual metrics, these gains are often accompanied by substantial deterioration elsewhere. For example, prioritizing low-confidence genes improves Pearson Δ and Cos. PCA but worsens Cos. LogFC and perturbation discriminability, suggesting that local uncertainty does not identify genes that provide useful context for subsequent generation. Our biological prior provides a strong directional generation scaffold. The biological-prior ordering introduced in D2R2D^2R^2 improves most metrics over Random, with only a modest decrease in Cos. LogFC, showing that regulatory structure alone provides an effective scaffold for gene generation. Further, reversing the prior order causes every metric to fall below Random, revealing that following regulatory dependencies in the wrong direction is actively detrimental rather than merely uninformative. Thus, the prior contributes not merely by imposing structure, but by organizing generation along a predictive regulatory direction. RPM turns static prior into perturbation-specific generation policies. Even with random initialization, RPM improves multiple metrics over Random, showing that response-level reinforcement learning can discover useful generation orders. With biological initialization, RPM converts the fixed prior into perturbation-specific policies that improve all three perturbation-effect metrics, with only modest trade-offs in centroid- and population-level agreement. Biological Analysis RPM preserved the generation prioritization of regulatory genes. As shown in Fig. 4a, regulatory genes were strongly enriched at early generation steps and declined thereafter, whereas the proportion of regulated genes progressively increased across steps. This result suggests that RPM tends to establish upstream regulatory context before generating downstream responses. RPM prioritizes perturbation-specific TFs to earlier generation steps. Next, we test whether RPM captures perturbation-specific regulatory rewiring. Perturbation-specific TFs are those inactive in unperturbed cells but acquire regulatory activity after perturbation (details in the Appendix). We quantify the RPM-induced change in generation priority as ΔOrder=OrderBio-prior−OrderRPM =Order_Bio -prior-Order_RPM, where a positive value indicates that RPM generates the gene earlier. In Norman19, perturbation-specific TFs exhibit significantly higher ΔOrder than non-specific TFs (Fig. 4b), demonstrating that RPM preferentially advances perturbation-specific regulators. RPM-promoted genes capture perturbation-specific transcription responses. Finally, we found that RPM-promoted genes (genes with ΔOrder>0 >0) were predominantly perturbation-responsive. As shown in Fig. 4c, RPM-promoted genes captured most perturbation-specific differentially expressed (DE) genes, and this pattern was consistent across all perturbations. RPM-promoted genes were also strongly enriched in cell-cycle-related pathways (see the Appendix), a major axis of perturbation response, especially given that many perturbations in Norman19 directly affect the proliferative state of cell. By prioritizing genes that characterize the perturbed cell state, RPM provides informative context for predicting the expression of later-generated genes. Figure 4: Biological analysis. (a) Proportions of regulatory and regulated genes across generation steps under the RPM-refined generation order. (b) Boxplot comparing ΔOrder between perturbation-specific TFs and non-specific TFs. Mann-Whitney U test, P=1.6×10−15P=1.6× 10^-15. (c) Boxplot showing the fraction of DE genes whose generation order was promoted by RPM or remained unpromoted. Each dot represents one perturbation, and paired fractions from the same perturbation are connected by lines. Paired Wilcoxon test, P=1.7×10−9P=1.7× 10^-9. Conclusion We introduced D2R2D^2R^2, a framework that makes gene generation order explicit in single-cell perturbation prediction. MDDM progressively generates masked expression tokens, while RPM initializes the ordering policy from biological prior and adapts it to each perturbation through reinforcement learning. Across Norman19 and H1, D2R2D^2R^2 ranks first across all metrics on Norman19 and leads in most metrics on H1. Controlled ablations show that biologically informed ordering is more effective than random and uncertainty-based strategies. Biological analyses further reveal that RPM preserves upstream regulatory context while prioritizing perturbation-specific regulators and responsive genes. References Adduri et al. (2025) Adduri, A. K.; Gautam, D.; Bevilacqua, B.; Imran, A.; Shah, R.; Naghipourfar, M.; Teyssier, N.; Ilango, R.; Nagaraj, S.; Dong, M.; et al. 2025. Predicting cellular responses to perturbation across diverse contexts with State. BioRxiv, 2025–06. Austin et al. (2021) Austin, J.; Johnson, D. D.; Ho, J.; Tarlow, D.; and Van Den Berg, R. 2021. Structured denoising diffusion models in discrete state-spaces. Advances in neural information processing systems, 34: 17981–17993. Bahdanau et al. (2016) Bahdanau, D.; Brakel, P.; Xu, K.; Goyal, A.; Lowe, R.; Pineau, J.; Courville, A.; and Bengio, Y. 2016. An actor-critic algorithm for sequence prediction. arXiv preprint arXiv:1607.07086. Bereket and Karaletsos (2023) Bereket, M.; and Karaletsos, T. 2023. Modelling cellular perturbations with the sparse additive mechanism shift variational autoencoder. Advances in Neural Information Processing Systems, 36: 1–12. Bunne et al. (2024) Bunne, C.; Roohani, Y.; Rosen, Y.; Gupta, A.; Zhang, X.; Roed, M.; Alexandrov, T.; AlQuraishi, M.; Brennan, P.; Burkhardt, D. B.; et al. 2024. How to build the virtual cell with artificial intelligence: Priorities and opportunities. Cell, 187(25): 7045–7063. Bunne et al. (2023) Bunne, C.; Stark, S. G.; Gut, G.; Del Castillo, J. S.; Levesque, M.; Lehmann, K.-V.; Pelkmans, L.; Krause, A.; and Rätsch, G. 2023. Learning single-cell perturbation responses using neural optimal transport. Nature methods, 20(11): 1759–1768. Chang et al. (2022) Chang, H.; Zhang, H.; Jiang, L.; Liu, C.; and Freeman, W. T. 2022. Maskgit: Masked generative image transformer. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, 11315–11325. Cui et al. (2024a) Cui, H.; Wang, C.; Maan, H.; Pang, K.; Luo, F.; Duan, N.; and Wang, B. 2024a. scGPT: toward building a foundation model for single-cell multi-omics using generative AI. Nature methods, 21(8): 1470–1480. Cui et al. (2024b) Cui, Z.; Xu, T.; Wang, J.; Liao, Y.; and Wang, Y. 2024b. Geneformer: Learned gene compression using transformer-based context modeling. In ICASSP 2024-2024 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 8035–8039. IEEE. Dixit et al. (2016) Dixit, A.; Parnas, O.; Li, B.; Chen, J.; Fulco, C. P.; Jerby-Arnon, L.; Marjanovic, N. D.; Dionne, D.; Burks, T.; Raychowdhury, R.; et al. 2016. Perturb-Seq: dissecting molecular circuits with scalable single-cell RNA profiling of pooled genetic screens. cell, 167(7): 1853–1866. Freimer et al. (2022) Freimer, J. W.; Shaked, O.; Naqvi, S.; Sinnott-Armstrong, N.; Kathiria, A.; Garrido, C. M.; Chen, A. F.; Cortez, J. T.; Greenleaf, W. J.; Pritchard, J. K.; et al. 2022. Systematic discovery and perturbation of regulatory genes in human T cells reveals the architecture of immune networks. Nature Genetics, 54(8): 1133–1144. He et al. (2026) He, S.; Zhu, Y.; Tavakol, D. N.; Ye, H.; Lao, Y.-H.; Zhu, Z.; Xu, C.; Chauhan, S.; Garty, G.; Tomer, R.; et al. 2026. Squidiff: predicting cellular development and responses to perturbations using a diffusion model. Nature methods, 23(1): 65–77. Klein et al. (2025) Klein, D.; Fleck, J. S.; Bobrovskiy, D.; Zimmermann, L.; Becker, S.; Palma, A.; Dony, L.; Tejada-Lapuerta, A.; Huguet, G.; Lin, H.-C.; et al. 2025. CellFlow enables generative single-cell phenotype modeling with flow matching. bioRxiv, 2025–04. Lotfollahi, Wolf, and Theis (2019) Lotfollahi, M.; Wolf, F. A.; and Theis, F. J. 2019. scGen predicts single-cell perturbation responses. Nature methods, 16(8): 715–721. Lutz et al. (2023) Lutz, I. D.; Wang, S.; Norn, C.; Courbet, A.; Borst, A. J.; Zhao, Y. T.; Dosey, A.; Cao, L.; Xu, J.; Leaf, E. M.; et al. 2023. Top-down design of protein architectures with reinforcement learning. Science, 380(6642): 266–273. Nie et al. (2025a) Nie, S.; Zhu, F.; Du, C.; Pang, T.; Liu, Q.; Zeng, G.; Lin, M.; and Li, C. 2025a. Scaling up Masked Diffusion Models on Text. In The Thirteenth International Conference on Learning Representations. Nie et al. (2025b) Nie, S.; Zhu, F.; You, Z.; Zhang, X.; Ou, J.; Hu, J.; Zhou, J.; Lin, Y.; Wen, J.-R.; and Li, C. 2025b. Large language diffusion models. arXiv preprint arXiv:2502.09992. Norman et al. (2019) Norman, T. M.; Horlbeck, M. A.; Replogle, J. M.; Ge, A. Y.; Xu, A.; Jost, M.; Gilbert, L. A.; and Weissman, J. S. 2019. Exploring genetic interaction manifolds constructed from rich single-cell phenotypes. Science, 365(6455): 786–793. Olivecrona et al. (2017) Olivecrona, M.; Blaschke, T.; Engkvist, O.; and Chen, H. 2017. Molecular de-novo design through deep reinforcement learning. Journal of cheminformatics, 9(1): 48. Piran et al. (2024) Piran, Z.; Cohen, N.; Hoshen, Y.; and Nitzan, M. 2024. Disentanglement of single-cell data with biolord. Nature Biotechnology, 42(11): 1678–1683. Qian, Dong, and Guo (2025) Qian, L.; Dong, Z.; and Guo, T. 2025. Grow AI virtual cells: three data pillars and closed-loop learning. Cell Research, 35(5): 319–321. Ranzato et al. (2015) Ranzato, M.; Chopra, S.; Auli, M.; and Zaremba, W. 2015. Sequence level training with recurrent neural networks. arXiv preprint arXiv:1511.06732. Roohani, Huang, and Leskovec (2024) Roohani, Y.; Huang, K.; and Leskovec, J. 2024. Predicting transcriptional outcomes of novel multigene perturbations with GEARS. Nature Biotechnology, 42(6): 927–935. Roohani et al. (2025) Roohani, Y. H.; Hua, T. J.; Tung, P.-Y.; Bounds, L. R.; Yu, F. B.; Dobin, A.; Teyssier, N.; Adduri, A.; Woodrow, A.; Plosky, B. S.; et al. 2025. Virtual Cell Challenge: Toward a Turing test for the virtual cell. Cell, 188(13): 3370–3374. Runge et al. (2018) Runge, F.; Stoll, D.; Falkner, S.; and Hutter, F. 2018. Learning to design RNA. arXiv preprint arXiv:1812.11951. Settles (2009) Settles, B. 2009. Active learning literature survey. Shao et al. (2024) Shao, Z.; Wang, P.; Zhu, Q.; Xu, R.; Song, J.; Bi, X.; Zhang, H.; Zhang, M.; Li, Y.; Wu, Y.; et al. 2024. Deepseekmath: Pushing the limits of mathematical reasoning in open language models. arXiv preprint arXiv:2402.03300. Shu et al. (2021) Shu, H.; Zhou, J.; Lian, Q.; Li, H.; Zhao, D.; Zeng, J.; and Ma, J. 2021. Modeling gene regulatory networks using neural network architectures. Nature Computational Science, 1(7): 491–501. Song et al. (2025) Song, B.; Liu, D.; Dai, W.; Mcmyn, N. F.; Wang, Q.; Yang, D.; Krejci, A.; Vasilyev, A.; Untermoser, N.; Loregger, A.; et al. 2025. Decoding heterogeneous single-cell perturbation responses. Nature cell biology, 27(3): 493–504. Tang et al. (2023) Tang, W.; Liu, R.; Wen, H.; Dai, X.; Ding, J.; Li, H.; Fan, W.; Xie, Y.; and Tang, J. 2023. A general single-cell analysis framework via conditional diffusion generative models. bioRxiv, 2023–10. Theodoris et al. (2023) Theodoris, C. V.; Xiao, L.; Chopra, A.; Chaffin, M. D.; Al Sayed, Z. R.; Hill, M. C.; Mantineo, H.; Brydon, E. M.; Zeng, Z.; Liu, X. S.; et al. 2023. Transfer learning enables predictions in network biology. Nature, 618(7965): 616–624. Vaswani et al. (2017) Vaswani, A.; Shazeer, N.; Parmar, N.; Uszkoreit, J.; Jones, L.; Gomez, A. N.; Kaiser, Ł.; and Polosukhin, I. 2017. Attention is all you need. Advances in neural information processing systems, 30. Wu et al. (2024) Wu, Y.; Wershof, E.; Schmon, S. M.; Nassar, M.; Osiński, B.; Eksi, R.; Yan, Z.; Stark, R.; Zhang, K.; and Graepel, T. 2024. Perturbench: Benchmarking machine learning models for cellular perturbation analysis. arXiv preprint arXiv:2408.10609. You et al. (2018) You, J.; Liu, B.; Ying, Z.; Pande, V.; and Leskovec, J. 2018. Graph convolutional policy network for goal-directed molecular graph generation. Advances in neural information processing systems, 31. Yu et al. (2026) Yu, C.; Wang, C.; Liao, B.; and Wu, T. 2026. scdfm: Distributional flow matching model for robust single-cell perturbation prediction. arXiv preprint arXiv:2602.07103. Yuan et al. (2021) Yuan, B.; Shen, C.; Luna, A.; Korkut, A.; Marks, D. S.; Ingraham, J.; and Sander, C. 2021. CellBox: interpretable machine learning for perturbation biology with application to the design of cancer combination therapy. Cell systems, 12(2): 128–140. Appendix Appendix A Background of Virtual Cell Cell, the fundamental functional unit of life, is a complex entity. The intricate cellular behaviors emerge from highly coordinated interactions among genes, proteins, metabolites, and signaling pathways across multiple biological scales (Yang et al. 2025). Understanding how these molecular systems collectively determine cellular responses to environmental stimuli and perturbations has long been a central goal in biology and medicine. However, the complexity, heterogeneity, and context dependency of cellular systems make comprehensive experimental characterization extremely challenging. To better understand and predict cellular behaviors, researchers have increasingly turned toward the construction of virtual cell models to simulate, predict, and steer cell behavior (Adduri et al. 2025; Bereket and Karaletsos 2023; Bunne et al. 2023; Cui et al. 2024a; He et al. 2026; Klein et al. 2025). The overarching vision of virtual cell modeling is to create predictive and mechanistic representations of cellular systems that can simulate how cells respond to perturbations, transition across developmental dynamics, and reorganize regulatory programs under varying biological conditions. Bunne et al. (Bunne et al. 2024) proposed the definition of Artificial Intelligence Virtual Cell (AIVC), which is a multi-scale foundation model that learns universal representations of biological entities across molecular, cellular, and tissue scales. Within this framework, cellular systems are represented through learned embeddings derived from multi-modal biological measurements, including genomic sequences, transcriptomics, proteomics, imaging, and spatial omics data. These representations can then be manipulated by neural network-based “virtual instruments” that simulate cellular transitions, perturbation responses, and developmental dynamics. In essence, the virtual cell paradigm seeks to transform cellular biology into a programmable computational system capable of both predictive simulation and mechanistic inference. The development of virtual cells has substantial scientific and clinical significance. Accurate virtual cell systems could enable large-scale in silico experimentation that would otherwise be prohibitively expensive or technically infeasible in physical laboratories. Furthermore, virtual cells may provide a unified computational framework for studying how molecular interactions collectively give rise to emergent cellular phenotypes across diverse biological contexts. Such systems have the potential to accelerate functional genomics, uncover regulatory mechanisms underlying development and disease, guide therapeutic target discovery, and improve precision medicine through patient-specific response prediction. The rapid development of high-throughput experimental techniques has significantly advanced progress in virtual cell modeling. CRISPR-based perturb-seq (Dixit et al. 2016), perturbation proteomics (Qian et al. 2024), and single-cell imaging (Chandrasekaran et al. 2024) now allow researchers to measure cellular states at unprecedented scale and resolution. In particular, datasets generated from genetic and chemical perturbation experiments (Heumos et al. 2025) provide paired observations of pre- and post-perturbation cellular states, offering valuable resources for learning cellular response dynamics. Driven by advances in experimental technologies, artificial intelligence has become a central methodology for virtual cell modeling. Given an initial cell state together with perturbation conditions, existing approaches aim to generate the corresponding post-perturbation cellular response. Current efforts mainly focus on perturbation-response prediction using transcriptomic measurements like scRNA-seq, since gene expression profiles provide scalable, information-rich representations of cellular states and are currently the most widely available modality for large-scale perturbation datasets. Existing approaches typically formulate virtual cell modeling as a conditional generative learning problem. Variational autoencoders (Bereket and Karaletsos 2023), transformer architectures (Adduri et al. 2025), graph neural networks (Roohani, Huang, and Leskovec 2024), and more recently diffusion-based generative models (He et al. 2026) are commonly employed to learn the mappings between perturbations and transcriptomic responses. These models aim to capture the high-dimensional structure of cellular states while generalizing across unseen perturbations, cell types, and biological conditions. Despite rapid progress, building reliable virtual cell systems remains an open challenge. Cellular systems exhibit complex nonlinear regulatory interactions, dynamic rewiring of gene regulatory networks, substantial biological heterogeneity, and strong context dependency across developmental and disease states (Yang et al. 2025). Moreover, current models often prioritize predictive accuracy while lacking mechanistic interpretability and biologically grounded generation processes. Future virtual cell systems must not only generalize across modalities and perturbations, but also support mechanistic reasoning, uncertainty estimation, and biologically faithful simulation of cellular dynamics (Wei et al. 2026). Developing deep learning architectures that can effectively integrate biological priors with generative modeling therefore remains a critical direction for future virtual cell research. Appendix B Experimental Details Datasets, Preprocessing, and Splits For both Norman19 and H1, we retain the top 1000 highly variable genes and keep their identities fixed throughout tokenization and generation. Binning statistics are fitted on training data and reused for validation and test examples. Control cells provide the reference expression used to form perturbation-induced changes. Norman19 contains 287 perturbation conditions, including 131 combinatorial conditions. We use the combination-prediction split with splitter seed 42; held-out combinations are reserved for evaluation. H1 uses the benchmark’s saved held-out-cell split. The dataset split is fixed across repeated runs, and the test split is used only for final evaluation. Evaluation Metrics Following the PerturBench protocol cited in the main paper, metrics are computed for each perturbation condition and then averaged across evaluation conditions. Let ^p y_p, py_p, and ctrly_ctrl denote the predicted, observed, and control mean-expression profiles for condition p. • Pearson Δ is the Pearson correlation between ^p−ctrl y_p-y_ctrl and p−ctrly_p-y_ctrl. • Cos. LogFC is the cosine similarity between predicted and observed log2 fold-change vectors relative to control, using a pseudocount of 0.1. • Cos. LogFC Rank ranks the matched observed perturbation against the other conditions by LogFC cosine similarity. Lower values indicate better perturbation-level discrimination. • Cos. PCA is the cosine similarity between predicted and observed condition centroids after projection into the reference PCA space. • Sym. KL is the symmetrized KL divergence between predicted and observed cell populations in the same reference PCA space. Lower values indicate closer distributional agreement. Baseline Protocol Comparison protocol. We designed the comparison so that methods differ in model design rather than in the data or evaluator. All methods use the same 1000-gene prediction space, perturbation and control annotations, training/validation/test partitions, and held-out test conditions. Predicted cells from every method are converted to the same AnnData schema and evaluated by the same PerturBench-compatible code, with the same reference controls, reference PCA space, condition aggregation, and metric definitions. No method-specific metric implementation or test-set post-processing is used. We preserve model-specific input representations and training procedures from the released implementations rather than forcing all methods into the discrete representation used by D2R2D^2R^2. For methods integrated into the benchmark (SAMS-VAE, Biolord, and Cell Flow), we use its unified data-loading and training interfaces. For standalone methods, including STATE, Squidiff, and scDFM, we adapt only the input/output interface required by the common split and evaluator. Model selection and early stopping use validation behavior only; test metrics are computed after model selection. Early stopping or a fixed training horizon follows the corresponding implementation, because an identical number of updates is not comparable across VAE, Transformer, diffusion, and flow-based objectives. Table 3: Model-specific configurations used for the baseline methods. Method Configuration SAMS-VAE Latent size 128; batch 256; dropout 0.4; learning rate 3.03×10−53.03× 10^-5 Biolord Latent size 512; batch 1000; dropout 0.4; learning rate 1.67×10−41.67× 10^-4 STATE Hidden size 768; eight Transformer layers; 12 attention heads; learning rate 10−410^-4 Cell Flow Batch 2000; hidden size 2048 with three layers; learning rate 10−510^-5 Squidiff Three layers; batch 64; learning rate 10−410^-4; EMA 0.9999 scDFM Hidden size 512; batch 32; learning rate 10−510^-5 D2R2D^2R^2 Configuration The denoising predictor uses 50 expression bins and a 12-layer Transformer with hidden dimension 768, 12 attention heads, and a SwiGLU intermediate dimension of 3072. We use RMSNorm and attention and path dropout of 0.1. Norman19 uses 64 control tokens and two perturbation tokens, whereas H1 uses one control token and four perturbation tokens. At inference time, generation proceeds for 20 steps and updates 50 genes per step, with an MDDM sampling temperature of 1.0. MDDM training. We first train the denoising predictor with the masked-token objective described in the main paper. The optimization settings are summarized in Table 4; the dataset-specific learning rate is selected on the validation split. GRPO training. After MDDM training, its parameters remain fixed while the RPM is optimized with the GRPO objective defined in the main paper. Each group contains four generation rollouts for the same example, and the resulting group-relative advantages are used for one clipped policy update. The complete optimization settings are given in Table 5. Table 4: MDDM training configuration. Dataset-specific values are reported as Norman19 / H1. Parameter Value Batch size 128 Optimizer AdamW Weight decay 10−210^-2 Learning rate (Norman19 / H1) 10−410^-4 / 5×10−45× 10^-4 Warmup fraction 0.01 Gradient clipping 1.0 Precision bfloat16 mixed EMA decay 0.998 Table 5: GRPO training configuration for RPM optimization. Parameter Value Batch size 96 Group size 4 Optimizer AdamW Learning rate 3×10−63× 10^-6 Weight decay 10−410^-4 Warmup fraction 0.01 Precision bfloat16 mixed Gradient clipping 1.0 GRPO clipping threshold 0.15 KL coefficient 0.05 Optimization iterations 1 Update cycles 50 Selection of bin representative values. We explored two strategies for mapping generated discrete tokens back to continuous expression values: assigning each sampled token the midpoint of its corresponding bin, and taking the posterior expectation over all bin midpoints. The expectation-based readout produced similar or slightly better performance on mean-profile metrics, but substantially degraded distribution-sensitive metrics such as MMD-PCA and Sym. KL, likely because averaging across bins reduces the variability of the generated cell population. We therefore use the bin-midpoint readout throughout our experiments. This mapping is deterministic and introduces no additional trainable parameters. Repeated Runs and Statistical Protocol The split is fixed across runs. The five reported seeds are 0,1,2,3,42\0,1,2,3,42\; they control training randomness, including initialization, data shuffling, and stochastic optimization, rather than dataset construction. Results in the main-results table are mean ± sample standard deviation over five independent runs. For D2R2D^2R^2, the MDDM checkpoint is fixed and the adaptive RPM branch is independently initialized and optimized for each seed. The generation-order ablation likewise fixes the MDDM and the 20-step generation budget so that only the selection policy changes. Lower-is-better metrics are sign-reversed when reporting paired improvements so that a positive difference always denotes improvement. These condition-level intervals are distinct from the training-seed standard deviations in the main table. Compute Environment Experiments were run under Ubuntu 22.04 with Python 3.10, CUDA 12.4, PyTorch 2.6, Lightning 2.6, and Scanpy 1.11 on NVIDIA A100 80GB PCIe GPUs. MDDM training and RPM optimization use bfloat16 mixed precision. Metric computation uses the PerturBench-compatible evaluation implementation released with the code. Appendix C Biological Analysis Analysis of regulatory-gene composition across generation steps. We inferred a cell-type-specific gene regulatory network (GRN) from unperturbed control cells using DeepSEM (Shu et al. 2021). We then assigned the 1,000 genes generated by D2R2D^2R^2 to two mutually exclusive categories: regulatory genes, defined as genes with at least one outgoing edge in the GRN, and the remaining genes are defined as regulated genes. This classification yielded 303 regulatory genes and 697 regulated genes. For each generation step, we calculated the proportion of generated genes in each category and plotted these proportions across steps. Comparison of ΔOrder between perturbation-specific and non-specific TFs. In this analysis, we identified perturbation-specific transcription factors (TFs) from DoRothEA regulon activity (Garcia-Alonso et al. 2019) using perturbation-versus-control expression changes. DoRothEA regulons were restricted to TF-target pairs represented among the genes modeled by D2R2D^2R^2. For each perturbation, we calculated the z-scored expression change of each gene and tested whether the changes across each TF’s target genes deviated significantly from the background distribution. The resulting P values were adjusted using the Benjamini–Hochberg procedure. TFs with an FDR<0.1FDR<0.1 were classified as perturbation-specific, whereas all remaining tested TFs were classified as non-specific. For each TF, we calculated ΔOrder as its biological-prior order minus its RPM-refined order. Thus, ΔOrder>0 >0 indicates that RPM moved the TF to an earlier generation step. We compared ΔOrder between the two groups using a one-sided Mann–Whitney U test. Figure 5: GSEA analysis using the ΔOrder values of all genes as the ranking metric. The GSEA analysis shows that RPM-promoted genes are enriched in cell cycle pathway (a) and G2-M transition pathway (b). Analysis of ΔOrder among differentially expressed genes. Genes with ΔOrder>0 >0 were defined as RPM-promoted genes, whereas those with ΔOrder≤0 ≤ 0 were defined as RPM-unpromoted genes. For each perturbation, we performed differential expression analysis against unperturbed control cells using scanpy.tl.rank_genes_groups. The top 100 genes with the largest absolute differential expression scores were defined as DE genes for that perturbation. Then, we calculated the fractions of DE genes classified as RPM-promoted and unpromoted, and compared the paired fractions across perturbations using a one-sided Wilcoxon signed-rank test. Gene set enrichment analysis (GSEA) of ΔOrder . Gene set enrichment analysis was performed on the averaged ΔOrder across perturbations using gseapy.prerank with Reactome Pathways 2024. A positive normalized enrichment score (NES) indicates enrichment among genes with higher ΔOrder , corresponding to genes promoted earlier by RPM. Appendix D The Results of GSEA Analysis As shown in Fig. 5, the RPM-promoted genes are strongly enriched in cell-cycle-related pathways, a major axis of perturbation response. References Adduri et al. (2025) Adduri, A. K.; Gautam, D.; Bevilacqua, B.; Imran, A.; Shah, R.; Naghipourfar, M.; Teyssier, N.; Ilango, R.; Nagaraj, S.; Dong, M.; et al. 2025. Predicting cellular responses to perturbation across diverse contexts with State. BioRxiv, 2025–06. Bereket and Karaletsos (2023) Bereket, M.; and Karaletsos, T. 2023. Modelling cellular perturbations with the sparse additive mechanism shift variational autoencoder. Advances in Neural Information Processing Systems, 36: 1–12. Bunne et al. (2024) Bunne, C.; Roohani, Y.; Rosen, Y.; Gupta, A.; Zhang, X.; Roed, M.; Alexandrov, T.; AlQuraishi, M.; Brennan, P.; Burkhardt, D. B.; et al. 2024. How to build the virtual cell with artificial intelligence: Priorities and opportunities. Cell, 187(25): 7045–7063. Bunne et al. (2023) Bunne, C.; Stark, S. G.; Gut, G.; Del Castillo, J. S.; Levesque, M.; Lehmann, K.-V.; Pelkmans, L.; Krause, A.; and Rätsch, G. 2023. Learning single-cell perturbation responses using neural optimal transport. Nature methods, 20(11): 1759–1768. Chandrasekaran et al. (2024) Chandrasekaran, S. N.; Cimini, B. A.; Goodale, A.; Miller, L.; Kost-Alimova, M.; Jamali, N.; Doench, J. G.; Fritchman, B.; Skepner, A.; Melanson, M.; et al. 2024. Three million images and morphological profiles of cells treated with matched chemical and genetic perturbations. Nature Methods, 21(6): 1114–1121. Cui et al. (2024a) Cui, H.; Wang, C.; Maan, H.; Pang, K.; Luo, F.; Duan, N.; and Wang, B. 2024a. scGPT: toward building a foundation model for single-cell multi-omics using generative AI. Nature methods, 21(8): 1470–1480. Dixit et al. (2016) Dixit, A.; Parnas, O.; Li, B.; Chen, J.; Fulco, C. P.; Jerby-Arnon, L.; Marjanovic, N. D.; Dionne, D.; Burks, T.; Raychowdhury, R.; et al. 2016. Perturb-Seq: dissecting molecular circuits with scalable single-cell RNA profiling of pooled genetic screens. cell, 167(7): 1853–1866. Garcia-Alonso et al. (2019) Garcia-Alonso, L.; Holland, C. H.; Ibrahim, M. M.; Turei, D.; and Saez-Rodriguez, J. 2019. Benchmark and integration of resources for the estimation of human transcription factor activities. Genome research, 29(8): 1363–1375. He et al. (2026) He, S.; Zhu, Y.; Tavakol, D. N.; Ye, H.; Lao, Y.-H.; Zhu, Z.; Xu, C.; Chauhan, S.; Garty, G.; Tomer, R.; et al. 2026. Squidiff: predicting cellular development and responses to perturbations using a diffusion model. Nature methods, 23(1): 65–77. Heumos et al. (2025) Heumos, L.; Ji, Y.; May, L.; Green, T. D.; Peidli, S.; Zhang, X.; Wu, X.; Ostner, J.; Schumacher, A.; Hrovatin, K.; et al. 2025. Pertpy: an end-to-end framework for perturbation analysis. Nature Methods, 1–10. Klein et al. (2025) Klein, D.; Fleck, J. S.; Bobrovskiy, D.; Zimmermann, L.; Becker, S.; Palma, A.; Dony, L.; Tejada-Lapuerta, A.; Huguet, G.; Lin, H.-C.; et al. 2025. CellFlow enables generative single-cell phenotype modeling with flow matching. bioRxiv, 2025–04. Qian et al. (2024) Qian, L.; Sun, R.; Aebersold, R.; Bühlmann, P.; Sander, C.; and Guo, T. 2024. AI-empowered perturbation proteomics for complex biological systems. Cell Genomics, 4(11). Roohani, Huang, and Leskovec (2024) Roohani, Y.; Huang, K.; and Leskovec, J. 2024. Predicting transcriptional outcomes of novel multigene perturbations with GEARS. Nature Biotechnology, 42(6): 927–935. Shu et al. (2021) Shu, H.; Zhou, J.; Lian, Q.; Li, H.; Zhao, D.; Zeng, J.; and Ma, J. 2021. Modeling gene regulatory networks using neural network architectures. Nature Computational Science, 1(7): 491–501. Wei et al. (2026) Wei, Z.; Wang, Y.; Gao, Y.; Wang, S.; Li, P.; Si, D.; Gao, Y.; Wu, S.; Li, D.; Dong, K.; et al. 2026. Benchmarking algorithms for generalizable single-cell perturbation response prediction. Nature Methods, 23(2): 451–464. Yang et al. (2025) Yang, T.; Wang, Y.-Y.; Ma, F.; Xu, B.-H.; and Qian, H.-L. 2025. Build the virtual cell with artificial intelligence: a perspective for cancer research. Military Medical Research, 12(1): 4.