Paper deep dive
Control-Anchored Residual Flow Matching Conditioned on Gene Geometry for Virtual Cell Perturbation Modeling
Quanquan Li, Yihe Chi, Liuyang Song, Hongbo Zhang, Jingyu Li, Xidong Xi, Conghua Wei, Yijie Sun, Yu Chen, Xin Liu, Qi Hu, Jing Ke, Guitao Cao
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 91%
Last extracted: 8/10/2026, 3:23:57 AM
Summary
The paper introduces GeneGeoFlow, a virtual cell modeling framework that predicts single-cell transcriptional responses to unseen genetic perturbations and drug combinations. It addresses the limitation of existing graph-based models that conflate stable gene relationships with intervention-specific response pathways. GeneGeoFlow employs a 'structure-dynamics separation' approach, deriving multi-scale spectral coordinates from Gene Ontology and control-derived coexpression networks to create perturbation-conditioned gene geometry. This geometry conditions a control-anchored residual flow matching model, coupled with condition-wise optimal transport and a Delta-correlation objective to align predicted and observed expression shifts. The model achieves high performance on Norman and ComboSciPlex benchmarks.
Entities (13)
Relation Signals (9)
GeneGeoFlow → achievesscoreon → Norman additive benchmark
confidence 95% · GeneGeoFlow achieves Pearson Delta scores of 0.8979 on the Norman additive benchmark
GeneGeoFlow → achievesscoreon → ComboSciPlex
confidence 95% · 0.9088 on five held-out drug combinations in the fixed ComboSciPlex test split.
GeneGeoFlow → uses → Gene Ontology
confidence 95% · GeneGeoFlow derives multi-scale spectral coordinates from Gene Ontology and control-derived coexpression networks.
GeneGeoFlow → uses → Control-derived coexpression networks
confidence 95% · GeneGeoFlow derives multi-scale spectral coordinates from Gene Ontology and control-derived coexpression networks.
GeneGeoFlow → employs → Optimal Transport
confidence 90% · Condition-wise optimal transport couples unpaired control and perturbed populations for training
GeneGeoFlow → employs → Delta-correlation objective
confidence 90% · a Delta-correlation objective aligns the predicted and observed condition-level expression-shift directions.
Spectral Scale Router → partof → GeneGeoFlow
confidence 90% · Perturbation-conditioned Spectral Scale Router and Graph Source Router then select the relevant spectral scales and graph sources for each gene
Graph Source Router → →
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:A central task in virtual cell modeling is predicting single-cell transcriptional responses to unseen genetic perturbations and drug combinations, and biological networks provide valuable priors on gene relationships. Existing graph-based models commonly use the same network to structure gene representations and mediate intergene interactions, thereby implicitly treating stable associations as perturbation-response pathways. Gene Ontology and control-derived coexpression networks encode relatively stable relationships rather than intervention-specific response directions or magnitudes. We therefore propose GeneGeoFlow, which conditions a control-anchored residual flow on gene-wise geometry derived from biological networks to learn intervention-specific transcriptional responses. GeneGeoFlow derives multi-scale spectral coordinates from Gene Ontology and control-derived coexpression networks. A perturbation-conditioned, gene-wise gating module selects relevant structural scales and network sources, yielding intervention-specific gene geometry. The resulting geometry conditions a control-anchored residual flow without explicitly propagating target-derived signals along the graph. Condition-wise optimal transport couples unpaired control and perturbed populations for training, while a Delta-correlation objective aligns the predicted and observed condition-level expression-shift directions. GeneGeoFlow achieves Pearson Delta scores of 0.8979 on the Norman additive benchmark and 0.9088 on five held-out drug combinations in the fixed ComboSciPlex test split. These results support perturbation-conditioned gene geometry as an effective structural prior for intervention-specific response prediction, without conflating stable gene relationships with response propagation.
Tags
Links
- Source: https://arxiv.org/abs/2608.06824v1
- Canonical: https://arxiv.org/abs/2608.06824v1
Trouble viewing inline? Open PDF directly →
Full Text
64,046 characters extracted from source content.
Expand or collapse full text
Control-Anchored Residual Flow Matching Conditioned on Gene Geometry for Virtual Cell Perturbation Modeling Quanquan Li1, Yihe Chi1, Liuyang Song2, Hongbo Zhang2, Jingyu Li3, Xidong Xi1, Conghua Wei1, Yijie Sun4, Yu Chen1, Xin Liu1, Qi Hu1, Jing Ke5, Guitao Cao1 Abstract A central task in virtual cell modeling is predicting single-cell transcriptional responses to unseen genetic perturbations and drug combinations, and biological networks provide valuable priors on gene relationships. Existing graph-based models commonly use the same network to structure gene representations and mediate intergene interactions, thereby implicitly treating stable associations as perturbation-response pathways. Gene Ontology and control-derived coexpression networks encode relatively stable relationships rather than intervention-specific response directions or magnitudes. We therefore propose GeneGeoFlow, which conditions a control-anchored residual flow on gene-wise geometry derived from biological networks to learn intervention-specific transcriptional responses. GeneGeoFlow derives multi-scale spectral coordinates from Gene Ontology and control-derived coexpression networks. A perturbation-conditioned, gene-wise gating module selects relevant structural scales and network sources, yielding intervention-specific gene geometry. The resulting geometry conditions a control-anchored residual flow without explicitly propagating target-derived signals along the graph. Condition-wise optimal transport couples unpaired control and perturbed populations for training, while a Delta-correlation objective aligns the predicted and observed condition-level expression-shift directions. GeneGeoFlow achieves Pearson Delta scores of 0.89790.8979 on the Norman additive benchmark and 0.90880.9088 on five held-out drug combinations in the fixed ComboSciPlex test split. These results support perturbation-conditioned gene geometry as an effective structural prior for intervention-specific response prediction, without conflating stable gene relationships with response propagation. 1 Introduction The space spanned by genetic perturbations and drug combinations far exceeds what single-cell experiments can measure (Norman et al. 2019; Lotfollahi et al. 2023). Moreover, destructive single-cell sequencing cannot provide paired pre- and post-intervention observations of the same cell, leaving available data sparse across conditions and unpaired across cells (Bunne et al. 2023). Virtual-cell models must therefore learn extrapolatable intervention responses from limited population-level observations. Consider the held-out Panobinostat+SRT1720 combination in Figure 1: starting from a control-cell population, the model must use response patterns learned from other training conditions to predict the combination-specific response distribution and gene-expression shift. Thus, the central problem is to learn from sparse, unpaired populations a response mapping that generalizes to unseen interventions. To enable such cross-intervention extrapolation, existing virtual-cell models learn transferable perturbation representations (Lotfollahi et al. 2019, 2023), transport unpaired cell populations (Bunne et al. 2023; Driessen et al. 2026), or model response distributions with flow, bridge, or diffusion dynamics (Klein et al. 2025; Chi et al. 2026; Yu et al. 2026; Yuan et al. 2026). Graph-based predictors such as GEARS further introduce biological priors, but use gene-relation graphs both to structure gene representations and to mediate intergene interactions (Roohani et al. 2024). Yet GO annotations encode curated functional relations (Ashburner et al. 2000), whereas control-derived coexpression reflects statistical associations among unperturbed cells; neither specifies the direction or magnitude of transcriptional change under a particular intervention. Using the same relations as response paths therefore conflates stable association with intervention-specific propagation. Thus, the key question is not how to propagate perturbations over a graph, but how to use biological networks as structural context while learning response dynamics from data. Answering this question requires what we term structure–dynamics separation: biological networks should provide gene-wise structural context, whereas intervention-specific response dynamics should be learned from data. This separation presents three technical challenges. First, a fixed graph source or spectral scale imposes the same structural view on every intervention, ignoring the condition-dependent roles of GO and control-derived coexpression. Second, control and perturbed cells have no one-to-one correspondence; arbitrary pairing confounds intervention effects with differences in baseline cell state (Bunne et al. 2023; Chi et al. 2026). Finally, cell-wise learning objectives can fit plausible individual cells without preserving the direction of condition-level gene-expression changes (Ahlmann-Eltze et al. 2025; Viñas Torné et al. 2026). The model must therefore jointly support intervention-adaptive structure selection, training couplings between unpaired populations, and condition-level response alignment. Guided by this principle, we introduce GeneGeoFlow, which operationalizes structure–dynamics separation by using biological networks as gene-wise structural coordinates and an independent response field to learn intervention-specific transcriptional changes. The model first extracts multi-scale spectral coordinates from GO and control-derived coexpression graphs. Perturbation-conditioned Spectral Scale Router and Graph Source Router then select the relevant spectral scales and graph sources for each gene, yielding gene geometry. Next, condition-wise optimal transport constructs training couplings between unpaired control and perturbed populations, and a control-anchored gene-wise residual flow learns intervention-specific population transitions conditioned on gene geometry. Finally, a Delta-correlation objective aligns the directions of predicted and observed condition-level expression changes. We evaluate GeneGeoFlow primarily on the Norman strict genetic-perturbation holdout protocol and use the fixed ComboSciPlex drug-combination split for cross-domain validation. On Norman strict holdout, gene geometry improves the five-fold mean Pearson Delta from 0.14500.1450 to 0.81530.8153 and outperforms the graph-free model in every fold; perturbation-conditioned geometry further improves over static geometry in matched folds, whereas explicit graph propagation yields no stable gain. On the established Norman additive benchmark, GeneGeoFlow achieves 0.89790.8979 Pearson Delta. On five held-out ComboSciPlex drug combinations, GeneGeoFlow achieves 0.90880.9088 Pearson Delta, compared with 0.86360.8636 for the explicit-propagation variant and 0.51690.5169 for the graph-free model. Within these matched comparisons, the dominant improvement arises from perturbation-conditioned gene geometry, while the evaluated propagation operator provides no consistent additional gain. Our contributions are summarized as follows: 1. We identify a fundamental role entanglement in existing graph-based perturbation models: biological networks are used both to represent gene organization and to mediate intervention-response propagation. We introduce a structure–dynamics separation perspective, in which biological networks define gene-wise structural context while intervention-specific response dynamics are learned from data. 2. We propose GeneGeoFlow, which constructs perturbation-conditioned gene geometry through multi-scale routing over Gene Ontology and control-derived coexpression graphs and uses this geometry to condition a control-anchored, gene-wise residual flow. Condition-wise optimal transport and a Delta-correlation objective address unpaired cell populations and condition-level response alignment, respectively. 3. Through matched comparisons with graph-free, static-geometry, and explicit-propagation controls on Norman and ComboSciPlex, we demonstrate that perturbation-conditioned gene geometry provides substantial improvements in predicting responses to unseen genetic perturbations and drug combinations. 2 Related Work Perturbation-conditioned gene geometry. Graph-based perturbation modeling has evolved from explicit relation propagation to module-aware and graph-conditioned response prediction. GEARS propagates perturbation information over gene relations (Roohani et al. 2024), while scBIG organizes gene responses through module-inductive representations (Ruan et al. 2026). scDFM further incorporates a graph-aware backbone into distributional flow matching (Yu et al. 2026). In parallel, scGen, CPA, chemCPA, biolord, AdaPert, and PerturbedVAE encode interventions through transferable, compositional, disentangled, or adaptive representations (Lotfollahi et al. 2019, 2023; Hetzel et al. 2022; Piran et al. 2024; Piao et al. 2026; Jiang et al. 2026). These methods primarily use graph-derived relations within response computation or adapt cell- and condition-level representations without modifying gene-wise structural coordinates. Spectral encoders such as SignNet provide global node coordinates but remain intervention agnostic (Lim et al. 2023). In contrast, GeneGeoFlow derives multi-scale spectral coordinates from GO (Ashburner et al. 2000) and control-derived coexpression, then uses perturbation-conditioned, gene-wise gates to select relevant scales and network sources. The resulting geometry provides intervention-specific structural context without propagating target signals. Geometry-conditioned residual flow. Unpaired population modeling has developed along two main lines. Optimal-transport methods construct population couplings without cell-level correspondences: CellOT and CMonge learn conditional control-to-treatment maps, while CINEMA-OT reduces treatment-associated variation before matching (Bunne et al. 2023; Driessen et al. 2026; Dong et al. 2023). Flow-based methods instead represent population change through a generative velocity field. CellFlow combines OT coupling with conditional flow matching, whereas scDFM introduces graph-aware structure into distributional flow matching (Klein et al. 2025; Yu et al. 2026). Despite their different parameterizations, these methods primarily address how a control population should be transported to a perturbed distribution, without explicitly separating stable gene structure from intervention-specific response dynamics. GeneGeoFlow restricts OT to endpoint construction and conditions a separate control-anchored, gene-wise residual field on ZpZ_p. Delta-correlation further aligns the predicted and observed condition-level expression shifts. Figure 1: Overview of GeneGeoFlow. (A) Condition-wise optimal transport couples unpaired control and perturbed populations to construct training endpoints. (B) Multi-scale spectral coordinates from the GO and control-derived coexpression graphs are selected by the Spectral Scale Router and fused by the perturbation-conditioned Graph Source Router, yielding gene-wise geometry ZpZ_p. (C) The control-anchored residual path interpolates between each coupled control–perturbed pair to form the intermediate state xtx_t. (D) A shared local projection maps (xt,xc,sp,Zp)(x_t,x^c,s_p,Z_p) to gene-wise features; the condition context modulates three FiLM residual blocks, and the velocity head outputs one scalar velocity per gene. (E) Integrating the learned cell-wide velocity field predicts the perturbed population from control cells. Panobinostat+SRT1720 is shown only as a running illustration: its 0.96610.9661 Pearson Delta is condition specific, whereas aggregate claims are evaluated across all held-out conditions. 3 Problem Formulation We define unseen-perturbation prediction as conditional population generation from unpaired single-cell observations. Let G be the number of modeled genes, xc∼ℙcx^c _c a control cell, and yp∼ℙpy^p _p a cell observed after perturbation p, all in ℝGR^G. Training data consist of finite control samples c=xici=1NcX_c=\x_i^c\_i=1^N_c and perturbed samples p=yjpj=1NpY_p=\y_j^p\_j=1^N_p for p∈trainp _train; destructive sequencing provides no cell-level correspondence between them. Each condition is described by a nonempty gene-aligned target mask sp∈0,1Gs_p∈\0,1\^G, with p=g:sp,g=1T_p=\g:s_p,g=1\. Let ℛ=GO,CER=\GO,CE\ and =W(r)r∈ℛG=\W^(r)\_r denote the fixed GO and control-coexpression graph priors. Training and test conditions satisfy train∩test=∅P_train _test= . For an unseen p⋆∈testp _test, the model receives only cX_c, sp⋆s_p , and G, with no sample from ℙp⋆P_p available for fitting. We seek a conditional generator y^p y^p ∼qθ(⋅∣xc,sp,), q_θ(\,· x^c,s_p,G\,), (1) ℙ^θ,p P_θ,p =∫qθ(⋅∣xc,sp,)dℙc(xc), = q_θ(\,· x^c,s_p,G\,)\,dP_c(x^c), where y^p y^p is one generated cell and ℙ^θ,p P_θ,p is the predicted population after marginalizing over control cells. For an unseen condition, the target is ℙ^θ,p⋆≈ℙp⋆ P_θ,p _p . Only the marginalized population ℙ^θ,p P_θ,p is identifiable from unpaired observations; the conditional coupling represented by qθq_θ is a modeling choice and is not interpreted as an individual counterfactual map. The model must therefore infer the response distribution of an unseen intervention from its target description and stable gene-structure priors, without assuming that those priors specify the response direction or magnitude. 4 Method GeneGeoFlow separates stable gene structure from intervention-specific response dynamics. It constructs perturbation-conditioned gene geometry ZpZ_p by routing multi-scale coordinates from fixed GO and control-derived coexpression graphs. Condition-wise OT supplies training endpoints, while a control-anchored gene-wise residual field uses ZpZ_p and perturbation context to predict responses; Delta correlation aligns condition-level shifts. The default model uses graph geometry only as structural context, whereas GeneGeoFlow-Prop adds target-signal propagation as a matched control. Figure 1 summarizes the workflow. 4.1 Perturbation-Conditioned Gene Geometry Fixed biological coordinates. We construct two weighted graphs over the modeled genes: a GO graph from biological-process annotation overlap (Ashburner et al. 2000) and a coexpression graph from absolute Pearson correlation across training control cells. Let gA_g be the GO term set of gene g and XcX^c the training-control expression matrix. Their edge weights are wghGO w_gh^GO =|g∩h||g∪h|, = |A_g _h||A_g _h|, (2) wghCE w_gh^CE =|Pearson(X:,gc,X:,hc)|, = |Pearson(X^c_:,g,X^c_:,h) |, where g and h index genes. Both graphs are sparsified, symmetrized, computed before training, and kept fixed; the exact sparsification rules are given in the supplementary material. Rather than treating raw graph edges as response pathways, we use spectral coordinates to summarize each gene’s higher-order structural position within each fixed graph. For source r∈ℛ=GO,CEr =\GO,CE\, let W(r)W^(r) be its adjacency and D(r)=diag(W(r))D^(r)=diag(W^(r)1). We obtain degree-normalized spectral coordinates from Lsym(r) L_sym^(r) =I−(D(r))−1/2W(r)(D(r))−1/2, =I-(D^(r))^-1/2W^(r)(D^(r))^-1/2, (3) Lsym(r)U(r) L_sym^(r)U^(r) =U(r)Λ(r), =U^(r) ^(r), Φ(r) ^(r) =(D(r))−1/2U(r), =(D^(r))^-1/2U^(r), where the columns of U(r)U^(r) are Laplacian eigenvectors, Λ(r) ^(r) contains their eigenvalues, and Φ(r) ^(r) gives the degree-normalized gene coordinates; we write λ(r)=diag(Λ(r))λ^(r)=diag( ^(r)) below. After discarding the first trivial mode, we retain the next 32 eigenpairs. We split them into low- and high-frequency blocks and apply separate SignNet-based Spectral Encoders (SPEs) (Lim et al. 2023), obtaining fixed gene-wise features blo,g(r)b_lo,g^(r) and bhi,g(r)b_hi,g^(r). The low-frequency block represents coarse graph organization, whereas the high-frequency block retains finer structural variation; keeping both provides candidate scales for perturbation-conditioned routing. Graph construction, isolated-node handling, and encoder dimensions are specified in the supplementary material. Perturbation-conditioned gene geometry. A fixed coordinate system cannot express that interventions may depend on different structural scales and graph sources. For the target set pT_p, let b¯g=[blo,g(r);bhi,g(r)]r∈ℛ b_g=[b_lo,g^(r);b_hi,g^(r)]_r and form the shared target-set embedding epgeo=|p|−1∑h∈pϕgeo(b¯h)e_p^geo=|T_p|^-1 _h _p _geo( b_h). The Spectral Scale Router (SSR) mixes the two frequency blocks within each source, and the Graph Source Router (GSR) then fuses the sources: αp,g(r) _p,g^(r) =sigmoid(frscale([blo,g(r);bhi,g(r);epgeo])), =sigmoid\! (f_r^scale([b_lo,g^(r);b_hi,g^(r);e_p^geo]) ), (4) zp,g(r) z_p,g^(r) =Mr(αp,g(r)blo,g(r)+(1−αp,g(r))bhi,g(r)), =M_r\! ( _p,g^(r)b_lo,g^(r)+(1- _p,g^(r))b_hi,g^(r) ), πp,g(r) _p,g^(r) =softmaxr(frsrc([zp,g(r);epgeo])), =softmax_r\! (f_r^src([z_p,g^(r);e_p^geo]) ), zp,g z_p,g =∑r∈ℛπp,g(r)zp,g(r), = _r _p,g^(r)z_p,g^(r), where ϕgeo _geo is a shared target-gene encoder, frscalef_r^scale and frsrcf_r^src are learned scalar scorers, and MrM_r maps each source to the common dzd_z-dimensional space. In Equation (4), αp,g(r) _p,g^(r) selects scale within each source, πp,g(r) _p,g^(r) fuses sources, and stacking the outputs gives Zp=[zp,1,…,zp,G]⊤∈ℝG×dzZ_p=[z_p,1,…,z_p,G] ^G× d_z. Thus, interventions share the same structural basis while selecting different structural scales and graph sources for each gene. The target mask conditions these selections and is not propagated along graph edges. 4.2 Geometry-Conditioned Gene-Wise Residual Flow Matching Given the perturbation-conditioned geometry ZpZ_p, the remaining task is to learn a conditional population transformation from unpaired control and perturbed observations. Condition-wise cell coupling and control-anchored path. For condition p, equal-size control and perturbed minibatches are coupled by discrete optimal transport (Dong et al. 2023; Klein et al. 2025; Chi et al. 2026): ρp⋆ _p =argminρ∈B∑i=1B‖xic−yρ(i)p‖22, = ρ∈ S_B _i=1^B \|x_i^c-y_ρ(i)^p \|_2^2, (5) yi y_i =yρp⋆(i)p, =y_ _p (i)^p, where B is the condition-specific minibatch size, B S_B is the set of index permutations, and ρp⋆ _p is the minimum-cost one-to-one assignment. Compared with arbitrary pairing, this assignment reduces baseline-state discrepancies between endpoints so that residual targets are less dominated by cross-cell mismatch. For each coupled pair, we perturb the measured control locally and construct a straight residual path: xi,0 x_i,0 =xic+σϵi, =x_i^c+σ _i, ϵi _i ∼(0,IG), (0,I_G), (6) xi,ti x_i,t_i =(1−ti)xi,0+tiyi, =(1-t_i)x_i,0+t_iy_i, ti t_i ∼(0,1), (0,1), ui u_i =yi−xi,0. =y_i-x_i,0. where σ is the fixed control-anchor noise scale, xi,tix_i,t_i is the sampled path state, and uiu_i is the constant target velocity along the straight path. Perturbing the control anchor encourages a locally smooth velocity field, while independent noise draws provide sample diversity when generating populations. The coupling and path are training constructions rather than lineage or individual counterfactual assignments. Condition context. The response backbone processes genes with shared parameters; within this backbone, the pooled Condition Context is the only channel through which information is shared across genes. It uses a separately parameterized perturbation embedding. It summarizes routed target coordinates as eprsp=|p|−1∑h∈pϕrsp(zp,h)e_p^rsp=|T_p|^-1 _h _p _rsp(z_p,h). A context encoder first constructs gene tokens and pools them into a shared condition representation: ηi,p,g _i,p,g =[xi,gc;zp,g;sp,g;eprsp], =[x_i,g^c;z_p,g;s_p,g;e_p^rsp], (7) τi,p,g _i,p,g =GeneTokenEncoder(ηi,p,g), =GeneTokenEncoder( _i,p,g), τ¯i,p τ_i,p =AttnPool(τi,p,gg=1G), =AttnPool(\ _i,p,g\_g=1^G), ci,p,ti c_i,p,t_i =[τ¯i,p;TimeEmbed(ti);eprsp], =[ τ_i,p;TimeEmbed(t_i);e_p^rsp], where ϕrsp _rsp is a response-specific shared target-gene encoder, GeneTokenEncoder maps each gene’s inputs to a token, AttnPool forms a cell-level summary, and TimeEmbed applies a sinusoidal encoding followed by an MLP. Learned-query attention pooling permits gene tokens to contribute unequally to the cell-level summary, while the time embedding identifies the current position along the control-to-perturbed path. The resulting ci,p,tic_i,p,t_i is the Condition Context in Figure 1; it is shared across genes within cell i but remains control-, perturbation-, and time-dependent. Local projection. For each gene, we concatenate the interpolated expression, control anchor, perturbation-conditioned geometry, and target state: ξi,p,ti,g _i,p,t_i,g =[xi,ti,g;xi,gc;zp,g;sp,g], =[x_i,t_i,g;x_i,g^c;z_p,g;s_p,g], (8) ℓi,p,ti,g _i,p,t_i,g =Ploc(ξi,p,ti,g)=SiLU(Wlocξi,p,ti,g+bloc). =P_loc( _i,p,t_i,g)=SiLU(W_loc _i,p,t_i,g+b_loc). The projection shares parameters across genes, while zp,gz_p,g retains gene-specific structural context. FiLM residual field and velocity head. The condition context first enters through an additive projection and then modulates residual blocks k=1,…,Kk=1,…,K, with K=3K=3. To simplify notation, write c=ci,p,tic=c_i,p,t_i and ℓg=ℓi,p,ti,g _g= _i,p,t_i,g, and suppress the common indices (i,p,ti)(i,p,t_i) inside the block recurrence. hg(0) h_g^(0) =ℓg+Pctx(c), = _g+P_ctx(c), (9) h~g(k) h_g^(k) =hg(k−1)+W2(k)SiLU(W1(k)LN(hg(k−1))), =h_g^(k-1)+W_2^(k)SiLU\! (W_1^(k)LN(h_g^(k-1)) ), [γ(k);β(k)] [γ^(k);β^(k)] =A(k)c+a(k), =A^(k)c+a^(k), hg(k) h_g^(k) =(1+γ(k))⊙h~g(k)+β(k), =(1+γ^(k)) h_g^(k)+β^(k), vθ,i,g v_θ,i,g =wout⊤hg(K)+bout, =w_out h_g^(K)+b_out, where PctxP_ctx, W1(k)W_1^(k), W2(k)W_2^(k), and A(k)A^(k) are learned maps; LNLN denotes layer normalization, ⊙ denotes element-wise multiplication, and γ(k)γ^(k) and β(k)β^(k) are condition-dependent FiLM scale and shift vectors. Here vθ,i,gv_θ,i,g abbreviates vθ,i,p,ti,gv_θ,i,p,t_i,g for the current (p,ti)(p,t_i), and stacking over genes gives vi=(vθ,i,g)g=1Gv_i=(v_θ,i,g)_g=1^G. Bias terms in the residual maps are suppressed for readability. FiLM treats c as a shared control signal that adapts gene-local feature processing at each block through condition-dependent scale and shift parameters. All residual-field and head parameters are shared over genes; genes interact only through the pooled Condition Context, not through graph propagation. Flow-matching and Delta-correlation objectives. The flow-matching term fits the OT-coupled residual velocity: ℒFM=1BG∑i=1B‖vi−ui‖22.L_FM= 1BG _i=1^B\|v_i-u_i\|_2^2. (10) This cell-level objective does not directly constrain the direction of the condition-level mean expression shift. We therefore use the noise-free endpoint residual to impose a Delta-correlation objective across genes: v¯p v_p =1B∑i=1Bvi, = 1B _i=1^Bv_i, d¯p d_p =1B∑i=1B(yi−xic), = 1B _i=1^B(y_i-x_i^c), (11) ℒΔ _ =1−corrκ(v¯p,d¯p). =1-corr_κ( v_p, d_p). where corrκcorr_κ is a Pearson correlation computed across the G genes with positive denominator floor κ. The complete objective is ℒ=ℒFM+λΔℒΔ,L=L_FM+ _ L_ , (12) where λΔ _ controls the strength of condition-level response alignment. The stabilized correlation operator and layer dimensions are given in the supplementary material. Explicit-propagation control. The default GeneGeoFlow uses graph spectra only as structural coordinates and never spreads the perturbation mask along graph edges. To isolate whether target-signal propagation adds value beyond this geometric conditioning, GeneGeoFlow-Prop introduces a matched learned low-rank spectral operator over the coexpression coordinates: rp(ℓ)=Φ(CE)diag(gψ,ℓ(λ(CE)))(Φ(CE))⊤sp.r_p^( )= ^(CE)diag\! (g_ψ, (λ^(CE)) )( ^(CE)) s_p. (13) where ℓ=1,…,C =1,…,C indexes the learned propagation channels, gψ,ℓg_ψ, is the filter for channel ℓ , λ(CE)λ^(CE) is the coexpression eigenvalue vector, and rp(ℓ)∈ℝGr_p^( ) ^G is the resulting gene-wise propagated channel. These channels are appended to the local input ξi,p,ti,g _i,p,t_i,g in Equation (8); all other components remain unchanged. This is a learned low-rank projection and reconstruction rather than an assumed causal diffusion process. Accordingly, the comparison evaluates this particular propagation operator under matched conditions rather than excluding all possible graph-propagation designs. Inference. Inference requires neither OT assignment nor perturbed cells. Given a control cell and an unseen target mask, the model constructs ZpZ_p and integrates ϵ ε ∼(0,IG), (0,I_G), (14) x(0) x(0) =xc+σϵ, =x^c+σε, dx(t)dt dx(t)dt =vθ(x(t),t,xc,sp,Zp),t∈[0,1]. =v_θ (x(t),t,x^c,s_p,Z_p ), t∈[0,1]. Integrating Equation (14) from t=0t=0 to t=1t=1 yields the terminal state x(1)=y^px(1)= y^p, one generated perturbed cell. Repeating this procedure over control cells and independent noise draws yields the predicted population ℙ^θ,p P_θ,p. Table 1: Norman additive results. Bold and underlined values denote the best and second-best results, respectively. Model MSE ↓ DE-Spearman ρ ↑ Pearson Δ ↑ Control 0.0184 N.A. N.A. Additive 0.0045 0.5564 0.9024 Geneformer 0.0041 0.3741 0.7732 GEARS 0.0139 0.5624 0.7421 CPA 0.0344 0.0713 0.3845 STATE 0.3006 0.5288 −0.0108-0.0108 CellFlow 0.0039 0.5503 0.8678 scDFM 0.0032 0.5705 0.8853 GeneGeoFlow 0.0030 0.5992 0.8979 5 Experiments 5.1 Experimental Setup We evaluated generalization to unseen interventions and used matched controls to isolate the effects of perturbation-conditioned gene geometry and explicit graph propagation. Datasets and evaluation protocols. Norman, our primary benchmark, contains K562 CRISPRa single- and two-gene interventions (Norman et al. 2019). Its five-fold strict-holdout protocol removes selected genes and every condition containing them from training, providing the main evidence for unseen genetic perturbations. We also report the established additive and single/double holdout protocols for comparison with scDFM. For cross-domain validation, we used ComboSciPlex, which contains A549 cells exposed to single agents and drug combinations (Lotfollahi et al. 2023). Its fixed split contains 24 training and seven held-out conditions; we focus on the five held-out drug–drug interventions and represent each combination by the union of its known targets. Complete preprocessing, split, and graph statistics are provided in the supplementary material. Baselines. Control reuses the control population and predicts no intervention shift. Additive sums available component-wise condition-mean residuals. Graph-free retains the same residual-flow backbone but removes spectral gene geometry, while GeneGeoFlow-Prop adds explicit target-signal propagation to GeneGeoFlow. The external tables additionally include scGPT (Cui et al. 2024), Geneformer (Theodoris et al. 2023), GEARS (Roohani et al. 2024), CPA (Lotfollahi et al. 2023), STATE (Adduri et al. 2025), CellFlow (Klein et al. 2025), and scDFM (Yu et al. 2026). For the Norman additive and single/double holdout comparisons, baseline values are transcribed from scDFM (Yu et al. 2026), while the GeneGeoFlow rows are obtained from our runs under the corresponding protocols. Because we did not rerun these external methods, the benchmark tables provide cross-source context, whereas matched local controls support mechanism attribution. Metrics and implementation. Our primary metrics were Pearson Delta, the gene-wise correlation between predicted and observed condition-mean residuals, and MSE, their numerical discrepancy. Population fidelity was additionally measured by RBF maximum mean discrepancy (MMD) and the Pearson correlation between predicted and observed gene-wise variance vectors. We averaged Norman results over data-split folds and ComboSciPlex results over conditions in its fixed test split unless stated otherwise. All locally trained models used the same preprocessing, evaluation interface, and training budget. Additional implementation and evaluation details are provided in the supplementary material. 5.2 Quantitative Results on Norman Figure 2: Matched controls isolate the contributions of biological gene coordinates and perturbation conditioning from explicit graph propagation. Panels (a)–(b) compare Graph-free, learned gene IDs, shuffled coordinates, and aligned target-set GeneGeoFlow on five Norman strict-holdout folds; in panel (a), hollow markers denote folds and filled diamonds denote five-fold means, while panel (b) reports mean MSE. Panel (c) reports Conditioned minus Static per Norman fold. Panel (d) reports GeneGeoFlow-Prop minus GeneGeoFlow for Norman folds and the ComboSciPlex mean. Each panel is a within-family comparison: panels (a)–(b) use target-set pooling, panel (c) uses the matched router run family, and panel (d) uses the dense full-mask response-encoder run family. Table 2: Norman single- and double-perturbation holdout results. Bold and underlined values denote the best and second-best results within each setting, respectively. Setting Model MSE ↓ DE-Spearman ρ ↑ Pearson Δ ↑ Single Control 0.0095 N.A. N.A. scGPT 0.0080 −0.1139-0.1139 0.4503 GEARS 0.0075 0.3569 0.6646 Geneformer 0.0036 0.3669 0.6955 CPA 0.0356 0.1168 0.2837 STATE 0.3333 0.6116 0.0004 CellFlow 0.0035 0.2860 0.7109 scDFM 0.0030 0.6957 0.7127 GeneGeoFlow 0.0029 0.7078 0.7421 Double Control 0.0207 N.A. N.A. scGPT 0.0153 −0.0665-0.0665 0.5693 GEARS 0.0156 0.2543 0.7552 Geneformer 0.0050 0.3468 0.7361 CPA 0.0357 0.3652 0.4176 STATE 0.3404 0.4071 0.0061 CellFlow 0.0049 0.5074 0.8095 scDFM 0.0047 0.5676 0.8357 GeneGeoFlow 0.0044 0.5918 0.8464 Table 2 provides benchmark context on the Norman single- and double-perturbation holdouts. GeneGeoFlow outperformed scDFM in both settings, reaching Pearson Delta scores of 0.74210.7421 and 0.84640.8464, with MSE values of 0.00290.0029 and 0.00440.0044, on the single- and double-perturbation holdouts, respectively. On the primary five-fold strict-holdout protocol, target-set GeneGeoFlow improved over Graph-free in every fold, increasing mean Pearson Delta from 0.14500.1450 to 0.81530.8153 and reducing MSE from 0.01670.0167 to 0.00390.0039. Figure 3: Qualitative population-level comparison on three fold-1 Norman strict-holdout conditions: ZBTB10+PTPN12 (a), BAK1+control (b), and DUSP9+PRTG (c). PCA and UMAP are fitted to a fixed stratified sample of real cells and then used to project both predicted populations. Points and contours denote cells and density; enlarged markers denote centroids. Aggregate population fidelity is evaluated separately in Table 4. The three conditions in Figure 3 span different levels of prediction difficulty, from partial recovery to close distributional agreement. Geometry-conditioned predictions follow the observed density more closely in panels (b) and (c), whereas panel (a) shows partial recovery; Table 4 provides the quantitative population-level test. Table 1 reports results under the established Norman additive protocol. GeneGeoFlow achieved 0.89790.8979 Pearson Delta, compared with 0.90240.9024 for Additive and 0.88530.8853 for scDFM, while also reporting the lowest MSE and the highest DE-Spearman ρ. 5.3 Quantitative Results on ComboSciPlex Table 3: ComboSciPlex results. Bold and underlined values denote the best and second-best results, respectively. Model MSE ↓ DE-Spearman ρ ↑ Pearson Δ ↑ Control 0.0324 N.A. N.A. scGPT 0.0031 −0.1261-0.1261 0.8322 CPA 0.0029 0.7906 0.8150 scDFM 0.0028 0.8289 0.8933 GeneGeoFlow 0.0021 0.8305 0.9088 Table 4: Population fidelity on held-out perturbations. Lower MMD and higher variance correlation are better; bold denotes the best result within each dataset. Dataset Variant MMD ↓ Var. Corr. ↑ Norman Graph-free 0.1488 0.3114 GeneGeoFlow-Prop 0.0965 0.6043 GeneGeoFlow 0.0821 0.6728 ComboSciPlex Graph-free 0.1639 0.2862 GeneGeoFlow-Prop 0.1011 0.5811 GeneGeoFlow 0.0894 0.6435 Table 3 places the ComboSciPlex result against reported baselines. GeneGeoFlow achieved the best results on all three reported metrics (MSE, DE-Spearman ρ, Pearson Δ ), including the strongest Pearson Delta (0.90880.9088 versus 0.89330.8933 for scDFM). In the matched local comparison, GeneGeoFlow improved Pearson Delta from 0.51690.5169 to 0.90880.9088 over Graph-free and reduced MSE from 0.07330.0733 to 0.00210.0021; explicit propagation reached 0.86360.8636 and 0.01980.0198. GeneGeoFlow improved all five unseen combinations by 0.19000.1900–0.48000.4800, including 0.27000.2700 to 0.71000.7100 on SRT2104+Alvespimycin. GeneGeoFlow also achieved the lowest MMD and the highest variance correlation on both datasets (Table 4). 5.4 Ablation Study Explicit-propagation control. Figure 2 separates structural coordinates from explicit target-signal propagation. Because the propagation analysis uses a dense full-mask response encoder, all conclusions in this panel are drawn within the same run family; its no-propagation reference therefore differs from the final target-set model. The evaluated low-rank operator produced mixed fold-wise effects on Norman, with a mean Pearson Delta change of +0.0101+0.0101, while decreasing ComboSciPlex performance by 0.04520.0452. It therefore provided no consistent benefit in this setting, without ruling out other propagation designs. Perturbation-conditioned geometry. Static geometry preserves the biological coordinates and target mask but removes the perturbation input from both routers. Conditioning improved every matched Norman fold, with a mean Pearson Delta gain of 0.03570.0357 (Figure 2(c)), isolating intervention-dependent scale and source selection. Gene coordinate alignment. With the response backbone fixed, Learned Gene ID, Shuffled Geometry, aligned geometry, and Graph-free reached mean Pearson Delta values of 0.39810.3981, 0.43870.4387, 0.81530.8153, and 0.14500.1450, respectively (Figure 2(a)). Their mean MSE values were 0.00910.0091, 0.00880.0088, 0.00390.0039, and 0.01670.0167 (Figure 2(b)), indicating that the gain cannot be explained by gene identity, input dimensionality, or arbitrary coordinate assignment alone. 6 Conclusion We introduced GeneGeoFlow to separate stable gene structure from intervention-specific response dynamics. It constructs perturbation-conditioned gene geometry from GO and control-derived coexpression networks to condition a control-anchored gene-wise residual flow; condition-wise OT and Delta correlation address unpaired populations and condition-level alignment. Matched comparisons on Norman and ComboSciPlex attributed the main gains to biologically aligned coordinates and perturbation-conditioned routing, whereas the evaluated low-rank propagation operator yielded no consistent benefit. These results support using biological networks as structural conditions for data-driven response fields. References A. K. Adduri, D. Gautam, B. Bevilacqua, A. Imran, R. Shah, M. Naghipourfar, N. Teyssier, R. Ilango, S. Nagaraj, M. Dong, C. Ricci-Tam, C. Carpenter, V. Subramanyam, A. Winters, S. Tirukkovular, J. Sullivan, B. S. Plosky, B. Eraslan, N. D. Youngblut, J. Leskovec, L. A. Gilbert, S. Konermann, P. D. Hsu, A. Dobin, D. P. Burke, H. Goodarzi, and Y. H. Roohani (2025) Predicting cellular responses to perturbation across diverse contexts with STATE. bioRxiv. External Links: Document, Link Cited by: §5.1. C. Ahlmann-Eltze, W. Huber, and S. Anders (2025) Deep-learning-based gene perturbation effect prediction does not yet outperform simple linear baselines. Nature Methods 22 (8), p. 1657–1661. External Links: Document, Link Cited by: §1. M. Ashburner, C. A. Ball, J. A. Blake, D. Botstein, H. Butler, J. M. Cherry, A. P. Davis, K. Dolinski, S. S. Dwight, J. T. Eppig, M. A. Harris, D. P. Hill, L. Issel-Tarver, A. Kasarskis, S. Lewis, J. C. Matese, J. E. Richardson, M. Ringwald, G. M. Rubin, and G. Sherlock (2000) Gene ontology: tool for the unification of biology. Nature Genetics 25 (1), p. 25–29. External Links: Document, Link Cited by: §1, §2, §4.1. C. Bunne, S. G. Stark, G. Gut, J. S. del Castillo, M. Levesque, K. Lehmann, L. Pelkmans, A. Krause, and G. Rätsch (2023) Learning single-cell perturbation responses using neural optimal transport. Nature Methods 20 (11), p. 1759–1768. External Links: Document, Link Cited by: §1, §1, §1, §2. C. Chi, Y. Huang, J. Xia, J. Zheng, Y. Liu, Z. Zang, and S. Z. Li (2026) Departures: distributional transport for single-cell perturbation prediction with neural Schrödinger bridges. Proceedings of the AAAI Conference on Artificial Intelligence 40 (25), p. 20544–20552. External Links: Document, Link Cited by: §1, §1, §4.2. H. Cui, C. Wang, H. Maan, K. Pang, F. Luo, N. Duan, and B. Wang (2024) scGPT: toward building a foundation model for single-cell multi-omics using generative AI. Nature Methods 21 (8), p. 1470–1480. External Links: Document, Link Cited by: §5.1. M. Dong, B. Wang, J. Wei, A. H. de O. Fonseca, C. J. Perry, A. Frey, F. Ouerghi, E. F. Foxman, J. J. Ishizuka, R. M. Dhodapkar, and D. van Dijk (2023) Causal identification of single-cell experimental perturbation effects with CINEMA-OT. Nature Methods 20 (11), p. 1769–1779. External Links: Document, Link Cited by: §2, §4.2. A. Driessen, D. A. Rajwade, B. Harsanyi, M. Rapsomaniki, and J. Born (2026) Conditional Monge gap enables generalizable single-cell perturbation modelling. Nature Machine Intelligence 8 (6), p. 984–996. External Links: Document, Link Cited by: §1, §2. L. Hetzel, S. Boehm, N. Kilbertus, S. Günnemann, M. Lotfollahi, and F. J. Theis (2022) Predicting cellular responses to novel drug perturbations at a single-cell resolution. In Advances in Neural Information Processing Systems, Vol. 35, p. 26711–26722. External Links: Link Cited by: §2. W. Jiang, Y. Liu, Y. Cai, E. Gao, J. Dong, E. Abbasnejad, L. Yao, and J. Q. Shi (2026) What makes a representation good for single-cell perturbation prediction?. In Forty-Third International Conference on Machine Learning, External Links: Link Cited by: §2. D. Klein, J. S. Fleck, D. Bobrovskiy, L. Zimmermann, S. Becker, A. Palma, L. Dony, A. Tejada-Lapuerta, G. Huguet, H. Lin, N. Azbukina, F. Sanchís-Calleja, T. Uscidda, A. Szalata, M. Gander, A. Regev, B. Treutlein, J. G. Camp, and F. J. Theis (2025) CellFlow enables generative single-cell phenotype modeling with flow matching. bioRxiv. External Links: Document, Link Cited by: §1, §2, §4.2, §5.1. D. Lim, J. D. Robinson, L. Zhao, T. E. Smidt, S. Sra, H. Maron, and S. Jegelka (2023) Sign and basis invariant networks for spectral graph representation learning. In The Eleventh International Conference on Learning Representations, External Links: Link Cited by: §2, §4.1. M. Lotfollahi, A. Klimovskaia Susmelj, C. De Donno, L. Hetzel, Y. Ji, I. L. Ibarra, S. R. Srivatsan, M. Naghipourfar, R. M. Daza, B. Martin, J. Shendure, J. L. McFaline-Figueroa, P. Boyeau, F. A. Wolf, N. Yakubova, S. Günnemann, C. Trapnell, D. Lopez-Paz, and F. J. Theis (2023) Predicting cellular responses to complex perturbations in high-throughput screens. Molecular Systems Biology 19 (6), p. e11517. External Links: Document, Link Cited by: §1, §1, §2, §5.1, §5.1. M. Lotfollahi, F. A. Wolf, and F. J. Theis (2019) scGen predicts single-cell perturbation responses. Nature Methods 16 (8), p. 715–721. External Links: Document, Link Cited by: §1, §2. T. M. Norman, M. A. Horlbeck, J. M. Replogle, A. Y. Ge, A. Xu, M. Jost, L. A. Gilbert, and J. S. Weissman (2019) Exploring genetic interaction manifolds constructed from rich single-cell phenotypes. Science 365 (6455), p. 786–793. External Links: Document, Link Cited by: §1, §5.1. Y. Piao, H. Kim, S. Kim, Y. Oh, J. Jeon, S. Hwang, J. Lim, W. Y. Kim, C. Park, and S. Ahn (2026) Learning adaptive perturbation-conditioned contexts for robust transcriptional response prediction. In Forty-Third International Conference on Machine Learning, External Links: Link Cited by: §2. Z. Piran, N. Cohen, Y. Hoshen, and M. Nitzan (2024) Disentanglement of single-cell data with biolord. Nature Biotechnology 42 (11), p. 1678–1683. External Links: Document, Link Cited by: §2. Y. Roohani, K. Huang, and J. Leskovec (2024) Predicting transcriptional outcomes of novel multigene perturbations with GEARS. Nature Biotechnology 42 (6), p. 927–935. External Links: Document, Link Cited by: §1, §2, §5.1. J. Ruan, R. Quan, L. Xu, Z. Yang, and Y. Yang (2026) Beyond independent genes: learning module-inductive representations for single-cell gene perturbation prediction. In Forty-Third International Conference on Machine Learning, External Links: Link Cited by: §2. C. V. Theodoris, L. Xiao, A. Chopra, M. D. Chaffin, Z. R. Al Sayed, M. C. Hill, H. Mantineo, E. M. Brydon, Z. Zeng, X. S. Liu, and P. T. Ellinor (2023) Transfer learning enables predictions in network biology. Nature 618 (7965), p. 616–624. External Links: Document, Link Cited by: §5.1. R. Viñas Torné, M. Wiatrak, Z. Piran, S. Fan, L. Jiang, S. A. Teichmann, M. Nitzan, and M. Brbić (2026) Systema: a framework for evaluating genetic perturbation response prediction beyond systematic variation. Nature Biotechnology 44 (6), p. 1050–1059. External Links: Document, Link Cited by: §1. C. Yu, C. Wang, B. Liao, and T. Wu (2026) scDFM: distributional flow matching model for robust single-cell perturbation prediction. In The Fourteenth International Conference on Learning Representations, External Links: Link Cited by: §1, §2, §2, §5.1. X. Yuan, X. Liu, Y. S. Zhang, Z. Zhang, H. Guo, and J. Tang (2026) PerturbDiff: functional diffusion for single-cell perturbation modeling. In Forty-Third International Conference on Machine Learning, External Links: Link Cited by: §1. Supplementary Material Appendix A Additional Method Details This supplement records the construction and implementation details omitted from the condensed method in the main paper. It uses the same notation as the main paper and introduces no additional model component. Graph Construction and Spectral Coordinates Fixed graph priors. Let the training control population be Xc=[x1c,…,xNcc]⊤∈ℝNc×GX^c=[x_1^c,…,x_N_c^c] ^N_c× G, and let gA_g be the set of GO biological-process terms assigned to gene g. The two graph weights are wghGO w_gh^GO =|g∩h||g∪h|, = |A_g _h||A_g _h|, (S1) wghCE w_gh^CE =|Pearson(X:,gc,X:,hc)|. = |Pearson(X^c_:,g,X^c_:,h) |. Equation (S1) is evaluated using training controls only. For GO, we retain the 20 highest-weight neighbors of each gene. For coexpression, we retain edges with absolute correlation at least 0.30.3 and augment them with the top 20 neighbors of each gene. The directed top-k matrices are symmetrized by arithmetic averaging, W←(W+W⊤)/2W←(W+W )/2, and fixed before model training. Perturbed cells are not used to construct either graph. GO annotations are read from the GAF 2.2 biological-process namespace using gene symbols; annotations with the NOT qualifier are excluded. The archived annotation file was generated on March 28, 2026 from the March 1, 2026 GO release. Spectral data. For r∈GO,CEr∈\GO,CE\, define D(r)=diag(W(r))D^(r)=diag(W^(r)1). The complete fixed spectral construction is Lsym(r) L_sym^(r) =I−(D(r))−1/2W(r)(D(r))−1/2, =I-(D^(r))^-1/2W^(r)(D^(r))^-1/2, (S2) Lsym(r)U(r) L_sym^(r)U^(r) =U(r)Λ(r),(U(r))⊤U(r)=I, =U^(r) ^(r), (U^(r)) U^(r)=I, Φ(r) ^(r) =(D(r))−1/2U(r), =(D^(r))^-1/2U^(r), =(Φ(r),λ(r))r∈GO,CE. =\( ^(r),λ^(r))\_r∈\GO,CE\. Equation (S2) is computed once for each graph source and cached. For isolated genes, the corresponding entry of (D(r))−1/2(D^(r))^-1/2 is set to zero. We discard the first trivial mode and retain the next 32 eigenpairs, ordered by increasing eigenvalue. The first 16 retained modes form the low-frequency block and the remaining 16 form the high-frequency block. Separate SignNet encoders resolve eigenvector sign ambiguity in the two blocks and produce blo,g(r)b_lo,g^(r) and bhi,g(r)b_hi,g^(r). The eigendecompositions are performed once and cached. Perturbation-Conditioned Routing Let the nonempty target set be p=g:sp,g=1T_p=\g:s_p,g=1\, and let b¯g=[blo,g(r);bhi,g(r)]r∈GO,CE b_g=[b_lo,g^(r);b_hi,g^(r)]_r∈\GO,CE\ collect the fixed multi-scale coordinates of gene g. The geometry embedding and source-specific routing variables are epgeo e_p^geo =1|p|∑h∈pϕgeo(b¯h), = 1|T_p| _h _p _geo( b_h), (S3) qp,g(r) q_p,g^(r) =[blo,g(r);bhi,g(r);epgeo], =[b_lo,g^(r);b_hi,g^(r);e_p^geo], αp,g(r) _p,g^(r) =sigmoid(frscale(qp,g(r))), =sigmoid (f_r^scale(q_p,g^(r)) ), z~p,g(r) z_p,g^(r) =αp,g(r)blo,g(r)+(1−αp,g(r))bhi,g(r), = _p,g^(r)b_lo,g^(r)+(1- _p,g^(r))b_hi,g^(r), zp,g(r) z_p,g^(r) =Mrz~p,g(r). =M_r z_p,g^(r). After aligning the source-specific spaces, the Graph Source Router computes ap,g(r) a_p,g^(r) =frsrc([zp,g(r);epgeo]), =f_r^src([z_p,g^(r);e_p^geo]), (S4) πp,g(r) _p,g^(r) =exp(ap,g(r))∑r′exp(ap,g(r′)), = (a_p,g^(r)) _r (a_p,g^(r )), zp,g z_p,g =∑rπp,g(r)zp,g(r), = _r _p,g^(r)z_p,g^(r), Zp Z_p =[zp,1,…,zp,G]⊤. =[z_p,1,…,z_p,G] . Equations (S3) and (S4) define the two-stage routing operation. The target set affects both gates through a shared coordinate encoder and mean pooling in epgeoe_p^geo. It does not alter the graph eigensystems or diffuse along graph edges. Residual Velocity Backbone Shared condition context. The geometry and response target-set encoders use separate parameters. For sample i, the response branch constructs eprsp e_p^rsp =1|p|∑h∈pϕrsp(zp,h), = 1|T_p| _h _p _rsp(z_p,h), (S5) ηi,p,g _i,p,g =[xi,gc;zp,g;sp,g;eprsp], =[x_i,g^c;z_p,g;s_p,g;e_p^rsp], τi,p,g _i,p,g =GeneTokenEncoder(ηi,p,g), =GeneTokenEncoder( _i,p,g), τ¯i,p τ_i,p =AttnPool(τi,p,gg=1G), =AttnPool(\ _i,p,g\_g=1^G), ci,p,ti c_i,p,t_i =[τ¯i,p;TimeEmbed(ti);eprsp]. =[ τ_i,p;TimeEmbed(t_i);e_p^rsp]. Equation (S5) defines the shared condition context. The GeneTokenEncoder contains two linear–normalization–GELU layers. Attention pooling supplies a shared cell-level context, while gene identity remains explicit in the local branch. Because ϕgeo _geo and ϕrsp _rsp are shared across genes, these embeddings remain defined for target genes that are absent as intervention targets during training. Local projection and FiLM residual field. For every gene, ℓi,p,ti,g _i,p,t_i,g =Ploc([xi,ti,g;xi,gc;zp,g;sp,g]), =P_loc([x_i,t_i,g;x_i,g^c;z_p,g;s_p,g]), (S6) vθ,i,g v_θ,i,g =VelocityHead(FiLMRes(ℓi,p,ti,g;ci,p,ti)), =VelocityHead (FiLMRes( _i,p,t_i,g;c_i,p,t_i) ), vi v_i =[vθ,i,1,…,vθ,i,G]⊤. =[v_θ,i,1,…,v_θ,i,G] . Here vθ,i,gv_θ,i,g abbreviates vθ,i,p,ti,gv_θ,i,p,t_i,g for the current condition and sampled time. The implemented backbone uses three residual blocks of width 256 and outputs one scalar velocity per gene. The default model contains no graph message-passing layer. Training Objective and Numerical Stabilization Within each condition, we form the squared-Euclidean cost matrix between sampled control and perturbed cells and solve an exact balanced linear assignment with the Hungarian algorithm. No entropic regularization or Sinkhorn approximation is used. For the resulting OT-coupled minibatch, the flow-matching term is ℒFM=1BG∑i=1B‖vi−ui‖22,ui=yi−(xic+σϵi).L_FM= 1BG _i=1^B\|v_i-u_i\|_2^2, u_i=y_i-(x_i^c+σ _i). (S7) Equation (S7) uses the noise-perturbed source only in the flow target. The Delta target excludes source noise: v¯p=1B∑i=1Bvi,d¯p=1B∑i=1B(yi−xic). v_p= 1B _i=1^Bv_i, d_p= 1B _i=1^B(y_i-x_i^c). (S8) Equation (S8) instead forms the condition-level target from noise-free residuals. For a,b∈ℝGa,b ^G, define (a) C(a) =a−1G(⊤a), =a- 1G(1 a)1, (S9) dκ(a,b) d_κ(a,b) =max(‖(a)‖2‖(b)‖2,κ), = (\| C(a)\|_2\| C(b)\|_2,κ ), corrκ(a,b) _κ(a,b) =⟨(a),(b)⟩dκ(a,b), = C(a), C(b) d_κ(a,b), ℒΔ _ =1−corrκ(v¯p,d¯p), =1-corr_κ( v_p, d_p), ℒ =ℒFM+λΔℒΔ. =L_FM+ _ L_ . Equation (S9) combines stabilized Delta correlation with flow matching. We use κ=10−8κ=10^-8; the denominator floor affects only numerically degenerate batches. Explicit-Propagation Control The default model uses the graph spectra only to form ZpZ_p. The matched GeneGeoFlow-Prop control additionally constructs C learned channels over the coexpression coordinates: rp(ℓ) r_p^( ) =Φ(CE)diag(gψ,ℓ(λ(CE))) = ^(CE)diag (g_ψ, (λ^(CE)) ) (S10) ⋅(Φ(CE))⊤sp,ℓ=1,…,C. ·( ^(CE)) s_p, =1,…,C. The channels in Equation (S10) are appended to the local features in Equation (S6); the geometry routers, context encoder, velocity backbone, training coupling, and objectives remain unchanged. Because Φ(CE) ^(CE) contains degree-normalized random-walk coordinates rather than an orthonormal Euclidean basis, this operator is a learned low-rank projection and reconstruction, not an exact causal diffusion process. Implementation and Evaluation Settings Datasets and preprocessing. We remove empty cells and genes, total-count normalize non-log-transformed matrices to 10410^4 counts per cell, and apply log1p. The modeled space contains the 5,000 highly variable genes selected with the Scanpy seurat procedure together with every perturbation target present in the dataset. Genetic interventions are represented by multi-hot target masks, while a drug combination uses the union of the mapped targets of its components. Table S1 reports the resulting modeled spaces and fixed graph sizes. Table S1: Dataset and fixed-graph statistics used in the reported runs. Edge counts refer to undirected edges after symmetrization. Dataset Genes Controls Cond. GO edges CE edges Norman 5,029 7,275 226 63,755 66,406 ComboSciPlex 5,028 1,451 31 30,365 50,818 Evaluation splits. Table S2 gives the condition counts used for optimization and evaluation. For the Norman strict holdout, each fold withholds 12 genes and assigns every condition containing at least one withheld gene to the test set; fold construction uses seed 42+fold index42+fold index. For ComboSciPlex, the released fixed split contains 24 training and seven test conditions, and one training condition is reserved for validation. The five reported drug–drug conditions are Panobinostat+Crizotinib, Panobinostat+Curcumin, Panobinostat+SRT1720, Panobinostat+Sorafenib, and SRT2104+Alvespimycin. The other two fixed-split test conditions are control+Alvespimycin and control+Dacinostat. Table S2: Numbers of training, validation, and test conditions after validation construction. The Norman additive counts apply to each of its five folds. Protocol Fold Train Val. Test Norman additive 1–5 180 9 37 Norman strict 1 172 9 45 2 180 10 36 3 185 10 31 4 179 9 38 5 166 9 51 ComboSciPlex fixed – 23 1 7 Table S3: Architecture settings used in the reported runs. Component Setting GO graph Biological process; top 20 Coexpression graph Control only; top 20; threshold 0.3 Retained spectra 32 per graph after first mode Frequency split 16 low / 16 high Target-set embeddings 32 / 32; mean pool; separate parameters Geometry width dzd_z 64 Gene-token width 128 Velocity backbone Width 256; 3 residual blocks Propagation channels 0 default; 8 in GeneGeoFlow-Prop Table S4: Optimization and sampling settings used in the reported runs. Component Setting Control noise σ 0.2 Optimizer AdamW; weight decay 10−510^-5 Learning rate 3×10−43× 10^-4 Schedule 2,000-step warmup; cosine decay Maximum steps 200,000 Per-condition sampling 256 cells Maximum batch size 512, grouped by condition Gradient clipping 1.0 Delta weight λΔ _ 0.03 EMA decay 0.999 Validation interval 5,000 steps Model selection Best validation Pearson Delta Inference solver Euler; 30 steps Benchmark export 128 controls; one draw per control Population evaluation 128 controls; 10 draws per control Hyperparameter exploration. On Norman fold 1 with seed 42, we evaluated λΔ∈0.01,0.03,0.05 _ ∈\0.01,0.03,0.05\, σ∈0.15,0.20,0.25σ∈\0.15,0.20,0.25\, propagation channel counts in 4,8,16\4,8,16\, propagation scales in 0,0.25,0.5,1.0\0,0.25,0.5,1.0\, and inference control counts in 128,256\128,256\. We selected the final settings by validation Pearson Delta and then fixed them for the reported runs. Reproducibility. We use seed 42 for Python, NumPy, and PyTorch, including condition splitting, control-cell sampling, and evaluation. The validation set is drawn from the training conditions only, and model selection uses the exponential-moving-average parameters at the checkpoint with the highest validation Pearson Delta. Because deterministic CUDA algorithms are not forced, exact bitwise equality can depend on the hardware and software stack. Computing environment. The primary experiments were executed on an Ubuntu 20.04.6 server with eight NVIDIA H20 GPUs (96 GB memory per GPU), two Intel Xeon Platinum 8480+ processors (224 logical CPU threads in total), and 1.8 TiB system memory. Each training process was bound to one H20 GPU, while independent folds and ablations were distributed across the eight devices. The software environment used Python 3.10.20, PyTorch 2.7.1 with CUDA 11.8 and cuDNN 9.1, NumPy 1.26.4, SciPy 1.15.3, AnnData 0.11.4, Scanpy 1.11.5, pandas 2.3.3, and scikit-learn 1.7.2; the server used NVIDIA driver 570.124.06. Evaluation scope and aggregation. The benchmark tables use the same 1,000-gene evaluation interface as their corresponding external protocols. Norman values are unweighted means over data-split folds, whereas ComboSciPlex values are fixed-split point estimates unless stated otherwise. Within an evaluation, condition-level metrics are first computed separately and then averaged with equal weight across conditions. Population MMD uses the biased estimator, including diagonal terms, with the summed RBF kernel k(x,y)=∑h∈ℋexp(−‖x−y‖22/(2h2))k(x,y)= _h (-\|x-y\|_2^2/(2h^2)) and ℋ=0.1,0.5,1.0,5.0H=\0.1,0.5,1.0,5.0\. Variance correlation compares predicted and observed gene-wise variances. Inference. OT assignments and perturbed cells are used only during training. At inference, each sampled control is initialized as x(0)=xc+σϵx^(0)=x^c+σε and advanced through the learned velocity field for 30 Euler steps; independent controls and noise draws form the generated population. Table S5: Comparison of reconstruction, differential-expression, and response-shift metrics on the Norman additive split. Bold and underlined values denote the best and second-best results, respectively. Model L2 ↓ MSE ↓ MAE ↓ DE-Spear. ρ ↑ Pearson Δ ↑ DS ↑ Pearson Δ ↑ Pearson Δ^20 _20 ↑ Control 3.9937 0.01839 0.03953 N.A. N.A. 0.5135 −0.1695-0.1695 −0.1297-0.1297 Additive 1.9395 0.00448 0.02276 0.5564 0.9024 0.9686 0.8584 0.9244 scGPT 3.4112 0.01349 0.03796 1.07×10−51.07×10^-5 0.5304 0.5404 0.2165 0.2414 Geneformer 1.9132 0.00410 0.02360 0.3741 0.7732 0.8241 −0.0078-0.0078 0.2239 GEARS 3.5531 0.01387 0.06624 0.5624 0.7421 0.8601 −0.0089-0.0089 0.2032 CPA 5.7629 0.03435 0.07894 0.0713 0.3845 0.6021 −0.0039-0.0039 0.2254 STATE 17.3330 0.30059 0.24705 0.5288 −0.0108-0.0108 0.5135 −0.0069-0.0069 0.2515 CellFlow 1.7064 0.00392 0.02207 0.5503 0.8678 0.9321 0.8395 0.8988 scDFM 1.7043 0.00315 0.02155 0.5705 0.8853 0.9737 0.8468 0.9260 GeneGeoFlow 1.6923 0.00303 0.02070 0.5992 0.8979 0.9775 0.8601 0.9260 Table S6: Comparison of reconstruction, differential-expression, and response-shift metrics on the Norman single- and double-perturbation holdout splits. Bold and underlined values denote the best and second-best results within each setting, respectively. Setting Model L2 ↓ MSE ↓ MAE ↓ DE-Spear. ρ ↑ Pearson Δ ↑ DS ↑ Pearson Δ ↑ Pearson Δ^20 _20 ↑ Single Control 2.6834 0.0095 0.0263 N.A. N.A. 0.5217 0.1618 0.1982 scGPT 2.5007 0.0080 0.0259 −0.1139-0.1139 0.4503 0.5680 0.0747 0.0798 GEARS 2.5641 0.0075 0.0466 0.3569 0.6646 0.8271 0.6356 0.7914 Geneformer 1.6962 0.0036 0.0191 0.3669 0.6955 0.8070 0.5620 0.6513 CPA 5.8060 0.0356 0.0853 0.1168 0.2837 0.5796 −0.0028-0.0028 0.0802 STATE 18.2543 0.3333 0.2693 0.6116 0.0004 0.5236 0.0154 0.2386 CellFlow 1.6758 0.0035 0.0191 0.2860 0.7109 0.8072 0.6138 0.6753 scDFM 1.6186 0.0030 0.0190 0.6957 0.7127 0.8914 0.6659 0.8116 GeneGeoFlow 1.5487 0.00291 0.01812 0.7078 0.7421 0.8956 0.6914 0.8197 Double Control 4.1882 0.0207 0.0423 N.A. N.A. 0.5322 −0.1303-0.1303 −0.0265-0.0265 scGPT 3.5171 0.0153 0.0362 −0.0665-0.0665 0.5693 0.5578 0.2814 0.2652 GEARS 3.7458 0.0156 0.0708 0.2543 0.7552 0.8766 0.6407 0.8413 Geneformer 2.0819 0.0050 0.0237 0.3468 0.7361 0.8067 0.6245 0.7261 CPA 5.7891 0.0357 0.0796 0.3652 0.4176 0.6311 0.2432 0.2870 STATE 18.4458 0.3404 0.2733 0.4071 0.0061 0.5289 −0.0023-0.0023 0.2580 CellFlow 2.1042 0.0049 0.0236 0.5074 0.8095 0.8622 0.6780 0.7155 scDFM 2.0309 0.0047 0.0235 0.5676 0.8357 0.9189 0.7769 0.8688 GeneGeoFlow 1.9732 0.00437 0.02306 0.5918 0.8464 0.9321 0.7867 0.8875 Table S7: Comparison of reconstruction, differential-expression, and response-shift metrics on the five held-out drug–drug conditions in the fixed ComboSciPlex split. Bold and underlined values denote the best and second-best available results, respectively; dashes indicate unavailable metrics. Model L2 ↓ MSE ↓ MAE ↓ DE-Spear. ρ ↑ Pearson Δ ↑ DS ↑ Control 5.3716 0.0324 0.0698 N.A. N.A. 0.5714 scGPT 1.6934 0.0031 0.0251 −0.1261-0.1261 0.8322 0.8571 CPA 1.6592 0.0029 0.0240 0.7906 0.8150 0.8980 scDFM 1.6567 0.0028 0.0220 0.8289 0.8933 0.8776 GeneGeoFlow 1.6491 0.0021 0.0213 0.8305 0.9088 0.8772 Baseline rows follow the same external benchmark sources as the condensed tables in the main paper, while the GeneGeoFlow rows report our runs under the corresponding protocols.