Generative Intervention Models for Causal Perturbation Modeling
Authors: Nora Schneider
Lars Lorch
Niki Kilbertus
Bernhard Schölkopf
Andreas Krause
Venue: ICML 2025 poster (2025)
Keywords: Causality
Causal Inference
Causal Modeling
Interventions
Perturbations
Generalization
Interpretability
Primary area: general_machine_learning->causality
TL;DR: We introduce a novel framework for predicting general perturbation effects  via causal models.

Abstract: We consider the problem of predicting perturbation effects via causal models. In many applications, it is a priori unknown which mechanisms of a system are modified by an external perturbation, even though the features of the perturbation are available. For example, in genomics, some properties of a drug may be known, but not their causal effects on the regulatory pathways of cells. We propose a generative intervention model (GIM) that learns to map these perturbation features to distributions over atomic interventions in a jointly-estimated causal model. Contrary to prior approaches, this enables us to predict the distribution shifts of unseen perturbation features while gaining insights about their mechanistic effects in the underlying data-generating process. On synthetic data and scRNA-seq drug perturbation data, GIMs achieve robust out-of-distribution predictions on par with unstructured approaches, while effectively inferring the underlying perturbation mechanisms, often better than other causal inference methods.

Generative Intervention Models for Causal Perturbation Modeling
Nora Schneider 12 Lars Lorch? 2
Niki Kilbertus
Bernhard Schölkopf 25 Andreas Krause?
Abstract
generalization to unseen instances requires an inductive bias (Ibragimov & Has' minskii, 1981). Perturbation effects are commonly predicted with either unstructured black-box models that predict the distribution shifts resulting from perturbations end-to-end or structured mechanistic models that characterize shifts as changes to an underlying causal data gettrratnnttiggg process.
We consider the problem of predicting perturba- tion effects via causal models. In many applica- tions, it is a priori unknown which mechanisms of system are modified by an external pertur- bation, even though the features of the pertur- bation are available. For example, in genomics, some properties of drug may be known, but not their causal effects on the regulatory pathways of cells. We propose generative intervention model (GIM) that learns to map these perturbation fea- tures to distributions over atomic interventions in a josttly-sstimaedd causal model. Contrary to prior approaches, this enables us to predict the distribution shifts of unseen perturbation features while gaining insights about their mechanistic ef- fects in the underlying gata-generating process. On synthetic data and scRNA-seq drug perturba- tion data, GIMs achieve robust out-of-distribution predictions on par with unstructured approaches. while effectively inferring the underlying pertur- bation mechanisms, often better than other causal inference methods.
Unstructured approaches commonly treat the task of pre- dicting perturbation effects as statistical learning problem, given the perturbed samples and the features of the perturba- tion (e.g., Hetzel et al., 2022; Lotfollahi et al., 2023; Bunne et al., 2023). For instance, when predicting drug effects, they use features like dosage levels or molecular footprints as covariates directly predict the perturbation response. This is practical because it allows making predictions for new perturbation features. However, these approaches do not ex- plicitly model the underlying data-generating mechanisms and how they change under a perturbation (Schölkopf, 2022), As a consequence, unstructured models are challenging to interpret, difficult to use in experiment design, and imprac- tical to use for perturbations with different feature spaces.
By contrast, mechanistic models characterize distribution shifts as interventions structured or causal model of the system (e.g., Pearl, 2009; Parascandolo et al., 2018; Gon- zalez et al., 2024; Roohani et al., 2024). This perspective provides explanations for the predicted shifts as it explic- itly infers the data-generatie mechanisms and how they change under an intervention. Crucially, a causal model generalizes new perturbations on the we know the atomic interventions on the system compo- nents induced by given perturbation. This assumption underlies most mechanistic approaches today and constrains them to settings where the perturbation targets are known (c.g., Maathuis et al., 2010) or, if unknown, inferred with- out considering the perturbation features (e.g., Mooij et al., 2020). By not leveraging the features of perturbation as covariates, causal models currently cannot make predictions for novel, non-atomic perturbations.
1. Introduction
Decision-making and experiment design in the sciences require predicting how systems respond to perturbations. Since testing all perturbations in a design space is often intractable, expensive, or unethical, scientists need data- driven approaches for predicting the effects of unseen per- turbations. For example, predicting how cells react to novel molecules can facilitate drug design in the vast space of candidate molecules (Sadybekov & Katritch, 2023; Tejada- Lapuerta et al., 2025). As in all learning tasks, however,
'Technical University of Munich and Helmholtz Munich, Ger- many 'Department of Computer Science, ETH Zurich, Switzer- land 'Munich Center for Machine Learning, Germany 'MPI for Intelligent Systems, Tübingen, Germany 'ELLIS Institute Tübingen, Germany. Correspondence to: Nora Schneider <nora.schneider @ helmholtz munich.de>.
In this work, we propose causal modeling approach for predicting the effects of general perturbations, whose fea- tures may not directly reveal how the perturbation modifies the causal generative process. To achieve this, we train a generative intervention model (GIM) that maps the pertur-
Proceedings of the 42nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Generative Intervention Models for Causal Perturbation Modeling
bation features to a distribution over atomic interventions in a jointly-lermed causal model. The intervened causal generative process then models the distribution after per- turbation. Since the GIM is shared across perturbations, it can systematically generalize to unseen perturbation fea- tures, similar to unstructured approaches, while providing mechanistic insights into how the data-generating process is altered. To evaluate the learned internal structure and per- turbation predictions separately, our experiments evaluate GIMs in causal discovery, intervention target identification, and end-to-end perturbation effect prediction on synthetic data. Especially in nonlinear systems, GIMs outperform most causal inference approaches in structure and interven- tion target discovery, while predicting distribution shifts on par with unstructured approaches in out-of-distribution settings. On single-cell RNA sequencing (scRNA-seq) data of drug perturbations (Srivatsan et al., 2020), we show that GIMs can generalize effectively to new drug dosages and characterize the distribution shifts on a selection of genes quantitatively and qualitatively more accurately than the baselines.
Lotfollahi et al., 2019), others generalizing to new pertur- bations by conditioning on y (Yu & Welch, 2022; Hetzel et al., 2022; Lotfollahi et al., 2023). Optimal transport meth- ods like CondOT (Bunne et al., 2022; 2023) learn transport maps and generalize by conditioning the map on 7 Domain- specific foundation models (Cui et al., 2024; Hao et al., 2024) learn similar conditioning tokens. These approaches provide little insight into the system and perturbation mech- anisms and are often highly specific to biochemistry.
Causal models When using causal modeling to predict perturbation effects, causal inference methods typically focus on learning a structured, causal model M from data collected under different perturbations D. Most methods assume that, for each perturbation, it is known how the system's 's causal mechanisms are modified. Under this assumption, M is estimated, e.g., through conditional independence tests (Spirtes et al., 2000) or model assump- tions (Chickering, 2002; Zheng et al., 2018; Yang et al., 2018: Lorch et al., 2021). However, these methods cannot be applied when the underlying atomic modifications are unknown. Recent approaches aim at additionally inferring the atomic modifications that underlie the observations D (Mooij et al., 2020; Brouillard et al., 2020; Hägele et al., 2023; Squires et al., 2020). However, these methods infer the changes only for the perturbations D and treat each perturbation independently, without leveraging the associated perturbation features 7. As result, they are unable to predict the effects of unseen perturbations y leaving p(x:"") undefined for general perturbations. To do so, existing methods would require knowing the specific atomic interventions in M caused by y"
2. Problem Statement and Related Work
We study system of variables x E IR with density p(x) under perturbations of its generative process. A perturba- tion of x is represented by a feature vector y E RR, which fully characterizes how p(x) shifts under the perturbation; we denote the resulting perturbed density as p(x: The perturbation features 7 are observable and can, for instance, represent the chemical properties of drug applied to cells or the context of an experiment. We observe data from K experimental environments, each defined by distinct per- turbation applied to the system. For each environment k, we observe the corresponding perturbation features y and the system's response, given by a data matrix (A) of i.i.d, sam- ples from p(x; e.g. the gene expression. Together, these observations form the full dataset D= = {(x (k) (*)) )
Hybrid mechanistic models Some works combine mech- anistic modeling with deep learning. However, these works either do not model perturbations (Parascandolo et al., 2018; Pervez et al., 2024) or are highly domain-specific, like Cel- 1Oracle (Kamimoto et al., 2023). GEARS (Roohani et al., 2024) or PDGrapher (Gonzalez et al., 2024), who rely on ex- pert knowledge to construct a graph neural network, which is not applicable in general settings.
Given our dataset and unseen perturbation y", aim to model
p(x|D;>)) D:y").
(L)
3. Generative Intervention Models
Our proposed framework leverages a causal model charac- terize the (posterior predictive) perturbed density p(x YD:)) for arbitrary 7 Before introducing our approach, we sum- marize relevant work on perturbation modeling and causal inference and discuss how they relate to our work.
In this section, we introduce generative intervention models, which form a novel class of densities p(x: 7) for perturba- tion modeling. Our approach leverages the inductive bias of causal model, whose data-generating process charac- terizes the density of under any perturbation y. To begin, we first describe causal models and how perturbations can be modeled as local modifications of the causal generative process. Then, we show how introducing an explicit genera- tive model of these local modifications allows us to link the perturbation features y to the causal model. This enables us to sample from p(x;"") for arbitrary perturbations y"
Autoencoders and optimal transport Common perturba- tion modeling approaches predict responses directly from y without a mechanistic data-generating process. Autoen- coders typically represent distribution shifts induced by per- turbations in a learned latent space, some only applying to combinations of known perturbations (Lopez et al., 2018;
2

---

Generative Intervention Models for Causal Perturbation Modeling
3.1. Causal Generative Process
particular, they independently infer separate I for each perturbation in D, without explicitly modeling how the fea- tures y relate to the intervention I. As result, they cannot infer the intervention I for new features and thus cannot make predictions for the corresponding system response in practical prediction tasks.
To model the densities y(x;")) we use a mechanistic model of the generative process of x. Specifically, we assume that x generated by a structural causal model (SCM; Pearl, 2009) with directed acyelic causal graph G E (0, } dxx and parameters 0 that encode the causal relationships among the variables x. The induced distribution p(x; G, 0) factorizes as product of independent causal mechanisms
Instead, we introduce an explicit generative process for the interventions I given the features 7. This generative inter- vention model (GIM) is fully parameterized by Φ and shared across perturbations. It characterizes the distributions of the intervention targets and parameters conditioned on y as
Ti Pa (r xG.; 0)
with p(x; G, = II= p(az Xa,; 0). Each conditional distribution models a variable o, given its direct causes (parents) XG,- Under this generative process, each 2; is in- dependent of its nondescendants in G given its parents XG.-
p(I:7, p(I: go (+)) (xb: (T,)))
(4)
where Is and he are arbitrary functions. Contrary to Hägele et al. (2023), we do not seek to directly infer I from D. In- stead, we view the interventions I as latent variables that contribute to the overall epistemic uncertainty given D and are marginalized out when modeling p(x Diy). Specif- ically, GIMs model the generative process of x under a perturbation y by marginalizing over its distribution of I as
Perturbations can be modeled as sparse mechanism shifts of the causal data-generating process, since perturbations typically only modify few causal mechanisms of a system at once (Peters et al., 2017; Schölkopf, 2022). To capture this, we model perturbations of x as local (atomic) interven- tions on a sparse subset of the conditionals in Equation (2). An intervention that targets a variable z, replaces its causal mechanism with an interventional mechanism
p(x; M, Φ) = MM)(I: dz
(5)
- /ʃ p(x vh: M) go (T))))))))))))))))) p ho (I, r)) dI dyb.
F2 (x; XG.; 0, x/)
with additional parameters vb. while not modifying any other conditionals of the factorization over G. We may model various types of interventions, such as hard or soft interven- tions, which differ in how they modify the observational conditionals (Appendix A). We denote the variable indices targeted by an intervention by E {0, 13dd and write
Unlike previous causal inference approaches (Mooij et al., 2020; Hägele et al., 2023), GIMs can thus generalize to sampling from p(x:7) for arbitrary unseen features 7.
3.3. Posterior Predictive Distribution of GIMs
Given a data collection D, our goal is to model the pre- dictive density p(x (D;7) for any (unseen) perturbation y (Section 2). From our generative model in Equation (5), we see that this requires posterior inference of the causal model M and GIM parameters Φ,, which map the perturbation fea- tures to the interventions in the causal model. For simplicity, we consider a crude approximation of the posterior, as indi- cated above, and simply learn maximum a posteriori (MAP) estimates of the causal model M and the GIM parameters Φ
M (G
I= (I.
For example, y could represent the molecular features of a drug and I its targets and effects in the generative process of the gene expressions x of cells, which characterized by M. In summary, p(x; M) denotes the causal generative model and p(x; MM, the model perturbed by the atomic intervention I. Next, we discuss how to bridge this approach to general perturbations y and ultimately model (Xx for arbitrary features
p(x Dii) = / p(x| M, Φ D) dM dæ
3.2. A Generative Process of the Interventions
p(x|MM :")
(6)
Perturbations of causal models can be modeled by atomic interventions Z. However, these atomic interventions are not observed in our problem setting, as we only observe the perturbation features 7. Thus, to model the general perturbations p(x;M, M)) the perturbation features y have to be linked to an atomic intervention I such that, intu- itively, p(x; M,T)= = p(x; M, for some I depending on y. Existing causal methods (c.g., Mooij et al., 2020; Hägele et al., 2023) directly infer I given samples from p)x;7). In
with MM,, 6" = arg maxM. p(M, Φ DD. We note that, since p(Z D.>) = I p(I:7, Φ) Do posterior in- ference for GIMs enables performing amortized posterior inference of I given D downstream (cf. Hägele et al., 2023). In the next section, we show how to infer the MAP estimates M", d= with gradient-based optimization.
Unlike unstructured, black-box models, our approach in- fers the causal graph and mechanisms M as well as the
3

---

Generative Intervention Models for Causal Perturbation Modeling
GIM parameters & approximate p(x D;7). As a re- sult, learning GIMs provides with mechanistic insights into the generative process of the system, allows general- izing to new perturbations and even enables performing multiple perturbations y jointly by generating their atomic interventions and applying them simultaneously M. In future work, GIMs may also quantify the uncertainty over its causal mechanisms (Lorch al., 2021), e.g., design the most informative next y: transfer its causal models to new applications that, e.g.. perturb M with known atomic interventions I directly; and incorporate priors about the causal structure G to improve its predictions (c.g., Roohani et al., 2024). None of these are straightforward black-box approaches (Section 2).
can express the density of observation under model M and an intervention I
d II po( Xc,;0)' pi( Xa, bb)"
p(x; M,
We describe two common choices for pú in Appendix A.
4.2. Maximum a Posteriori Estimation
With all causal model conditionals defined, we compute the likelihood of the data p(x; M,, 6) in Equation (5), Given the collection D of K datasets (X ,y()),))))) our goal described in Section 3.3 is to infer the MAP estimate
M " = arg max M.
D)
4. Inference
For this, we define suitable priors over M and & and then maximize the joint distribution over D. We model M and & as latent quantities, shared across the environments D, and the environments as conditionally independent given MM and & (Figure 1). Thus, the complete-data joint density is
In this section, we describe how to jointly infer the MAP estimates of the causal model M and the GIM parameters & with gradient-based optimization.
4.1. Model Components
Throughout this work, we make set of modeling choices to render the joint MAP inference of . and & tractable. We consider, e.g., linear Gaussian causal mechanisms given by
log log D) logr(M) log p(()
K
(9)
log logp(X (k)
JM, (()))
Pr xa,; 0) = (x;;; O,T (G, o x), exp(o:)) (7)
Here, the log likelihood for X (k) is obtained by summing the individual log likelihoods in Equation (5), since the samples are i.i.d. Before describing our priors and approximation of the likelihood, we first explain how we treat the discrete graph in M during estimation.
(e.g., Zheng et al., 2018) and nonlinear mechanisms with pr( xx,, (x;;;; MLPP, (G, o x), exp((()) (8)
(e.g., Brouillard et al., 2020; Lorch et al., 2021), where denotes elementwise multiplication and MLPe, denotes multi-layer perceptron (MLP) mapping Rd - IR with trainable parameters 0,. We learn the noise scale for each variable jointly to mitigate data scale artifacts (Reisach et al., 2021; Ormanicc et al., 2025). This then fully specifies the causal generative model p( (x; MM)
Differentiable inference of the causal model M We cannot perform gradient ascent on Equation (9) naively, because the graph matrix G in M is discrete and needs be acyclic (Section 3). To address this, representation of M by Lorch et al. (2021). Specifically, we use a latent representation Z E 2x2xxxz to model G probabilistically. where pz is a parameter controlling the rank of the graphs that the distribution can represent. The graph distribution defined as
the continuous
Similar to Brouillard al. (2020) and Hägele et al. (2023), model the targets I independent Bernouliis and the parameters ab as independent Gaussian distributions
p( (()))())) ( pl (bb ho ()) (n ,),1))
- 5 p( (Z)p(G Z)
p(G) -
d with PG = II II Bern ( ==1
The p.m.f. of Bern and p.d.f. of N are applied elementwise to each output of Jø and ha- The standard deviation 7)a is hyperparameter, and o-)) is the sigmoid function. In practice, 26 and hss are MLPs. Together with p(x; we
where > is an inverse temperature parameter. Here, and Zıy denote the i-th and j-th rows of the first and second (d x pz) matrices in Z, respectively. For pz 2 d, the con- ditional distribution p( G Z) can represent any adjacency matrix without self-loops. We denote this continuous model
'The specific choices in (7, 8) are not essential, long the dependence on G allows for Gumbel-softmax relaxation. In our scRNA-seg experiments, we use zero-inflated log.V / noise models.

---

Generative Intervention Models for Causal Perturbation Modeling
All Gaussian densities are again applied elementwise, and cGG) is the NO-BEARS acyclicity constraint (Lee et al., 2019), which is defined in Appendix D.4. The priors contain several hyperparameters concerning the sparsity (3z, 3M), regularization (πg 1/z, Tj. Aa), which require tuning on held-out data but can often be shared.
Perturlation Features
Intervestion Targets& & Parameters
System Response
Causal Model (Grrpph & Parzmeters)
(K)
I
X
M
(drug identity, dosage)
Igene
K
Stochastic optimization Finally, to optimize the joint p(M Φ D) with respect to Ma and Φ, we approx- imate the integrals left implicit in the likelihood and pri- ors. For this, we draw samples from both the graph model p(G | Ma) and the GIM p(Z: Φ) to compute the log like- lihood (10) and the log model prior (Equation (12)) via Monte Carlo integration. Overall, our MAP objective, the approximate continuous log joint density, is given by
(K experiments)
Φ
GIM Parameters tacural network weights)
Figure 1. Graphical model for the MAP estimation of the causal model M and GIM parameters Φb The atomic interventions I are marginalized out during inference. The dataset D consists of K pairs of observed perturbation features 7 7(h) and data matrices (()), For clarity, we depict (*) as observed random variable, but there is no prior over and we treat the perturbation features
log Dog((Mas $ D)
~ logp(Ma) + +lgp(()
K
constant throughout. Gray labels indicate concrete examples of the variables in the context of drug perturbation experiments.
+ +log
1
Σ ===1=
II p( (()), 0, 0, (m)(k)
M
with G(m) ~ p and [(m,k) p(I(k) (+) for m E {I....,M} The log likelihood estimator is biased due to the logarithm but consistent M - co. To compute p(Mo), we use the same approximation with independent samples from p(G G Ma). Since our causal generative process in Section 3.1 is only well-defined if G is acyclic, we follow Zheng et al. (2018) and maxi- mize log p(Ma. Φb D) with respect to M. and & with the augmented Lagrangian method, which iteratively updates the acyclicity penalties À and μ (Appendix D.4). When D contains observational data with p(x;7) = p(x), we set I = O, which implies M(x; = p(x: M). The MAP estimate for M is obtained via thresholding M. as G = ZZ Similarly, we use arg max, p 7 Φ) to sample from p(x "M"" at test time.
representation by Ma = {Z, 0). Under this generative model, the likelihood of the data D is given by
p(D Ma. Φ)
(10)
= -/NGIIM)
I( (8)
G, 0, Φ: ()) aG
(see Appendix B. The likelihood is the product of the GIM densities in Equation (5) over all observed perturba- tions in expectation over the graph given Ma. In this gen- erative model, expectations over M(Mas, Φ 'D) converge to those over p(M, & ()) with G = 1 1Z Z as a + oo (Lorch et al., 2021). Hence, instead of log (M, Φ. D), we optimize log P(Ma, Φ D) with respect to Ma.
Priors Our priors for learning GLMs capture the assump- tion that causal dependencies and interventions are sparse (Schölkopf, 2022). In addition to L2 and sparsity regular- ization factors in all priors, the priors penalize cycles in G. Specifically, we use the GIM parameter prior
To compute the model gradients through the Monte Carlo samples during optimization, the Gaussians modeling @ and s use the default reparameterizzation trick. The Bernoullis modeling I and G are reparameterized with the Gumbel- sigmoid (Maddison et al., 2017; Jang et al., 2017), except for the gradients with respect to 0, where we use the discrete G. In Appendix B, we provide the explicit forms of the gradients of log p(Ma. sb D) with respect to Ma and Φb
K
p(Φ) x
II exp ( E p (00
[ [l)))
targel sparsity
(11)
Identifiability Under standard assumptions on the func- tional form and noise models, identifiability results for causal discovery with unknown interventions (Brouillard et al., 2020) apply directly to GIMs. Specifically, in the large-sample limit, the MAP estimates of GIMs recover the true intervention targets in the training data and causal graph from the same interventional Markov equivalence class as the true graph (Yang et al., 2018). This result holds because, in the large-sample limit, the posterior is domi- nated by the likelihood, so the MAP estimates converge to
which encourages sparse interventions via an In penalty on the predicted targets. Additionally, we use the model prior
p(Ma) x NV 0, (O: 0, n8)
No; O, n
-exp
EppG 2) 3M G II.)))
mochanism sparsity
(12)
exp
Ep(/| 2)) (A G) + (c))))
acyclicity
5

---

Generative Intervention Models for Causal Perturbation Modeling
the same optimizers characterized in Brouillard et al. (2020). This further requires that the true interventional equivalence class contains graph in the support of the prior over causal graphs. Moreover, if mapping from perturbation features to intervention targets exists and is expressible by Φb and the prior over & supports this mapping, the true training targets can be recovered through optimization of Φ. Identifiability for unseen perturbations depends on the informativeness of the features 7, which we evaluate empirically Section 6.
e E {0.25, 0.75, 1.25, 1.75, 2.25) of the same intervention targets and Hill function mechanisms as in the training data. In the fully OOD task, we sample 20 new target and Hill function pairs distinct from the training set and perform perturbations at the training dosages e E (0.5, 1, 1.5,
We create the features Y by standardizing the concatenation of I and v. applying principal component analysis (PCA) across all training environments of system, and keeping the top 15 principal components. This representation later also allows us to study the effect of the information content in y on generalization. The data for a given y is sampled from the interventional distribution of the SCM p(x G, II tb).
5. Experimental Setup
We evaluate GIMs at predicting perturbation effects using both synthetic and real-world drug perturbation data.? GIMs estimate causal model . of x and jointly learn to predict which interventions I in M characterize p(x y): Al- though accurate estimates of these latent components are closely linked to the accuracy of p(x 7D; their full identi- fication is not strictly necessary for generalization. Thus, we first separately evaluate the inferred causal mechanisms and the predicted distributions on synthetic data with known gen- erative processes. On scRNA-seq drug perturbations, where these mechanisms are unknown, we focus on evaluating the predicted distribution shifts.
SciPlex3 drug perturbation data We also evaluate the predictive performance of GIMs on seRNA-seq data by Sri- vatsan al. (2020). This dataset measures gene expressions under drug perturbations at four different dosages (10n.M, 100nM, 1,M, and 10μM) across three human cancer cell lines (A549, K562, MCF7). For our experiments, we fo- cus on four drugs (Belinostat, Dacinostat, Givinostat, and Quinostat) from the epigenetic regulation pathway that were reported as among the most effective in this dataset and used in previous studies (Srivatsan al., 2020; Hetzel et al., 2022). We apply standard single-cell preprocessing steps detailed in Appendix D.2. Given the high dimensionality of the data, we focus our analysis on the top 50 marker genes, i.e., genes with strongest post-perture post-perturbation effects overall. We construct the perturbation features y by combining the one-hot encodings of the drugs and dosages as well as the dosage values. For evaluation, we create test sets by holding out, one at time, the highest dosage (10μM) of each drug, resulting in four unique training-test splits per cell type.
5.1. Datasets
Synthetic data To generate ground-truth systems, we cre- ate SCMs with additive Gaussian noise (Pearl, 2009). For this, we follow previous work (e.g.. Brouillard et al., 2020) and sample DAGs from either an Erdös-Rény: (ER) (Erdös et al., 1960) or a scale-free (SF) (Barabási & Albert, 1999) distribution with d = 20 nodes and 2d edges in expectation. We then generate random linear or nonlinear mechanisms. Appendix D.1 describes the synthetic data in detail.
5.2. Baselines
Generating synthetic perturbation mechanisms with features y that allow for generalization requires care. We draw in- spiration from biological simulators (Dibacinia & Sinha, 2020) and the biochemical perturbations evaluated later on and utilize et al., 2012) to model different perturbations of the ground- truth SCM. Hill functions are nonnegative, saturating func- tions that satisfy h(0)- = 0 and are used, e.g., to model transcription factor binding in gene regulation (Chu et al., 2009, see also Figure 2C). Inspired by drug compounds, each perturbation on the SCM consists of one or two ran- fixed but randomly-
We consider GIMs with different causal mechanisms. On the linear and nonlinear SCM data, GIMs learn Gaussian linear and MLP mechanisms as in Equations (7) and (8), respec- tively. On the ScRNA-seg count data, we use zero-inflated log- Dog-normal MLP mechanisms (see Appendix D.4). We com- pare GIMs with black-box and causal inference approaches (Section 2). Our baselines consist of causal discovery meth- ods for unknown interventions: BaCaDI (Hägele et al., 2023), GnIES (Gamella et al., 2022), UT-IGSP (Squires et al., 2020), and JCI-PC (Mooij et al., 2020; Spirtes et al., 2000). BaCaDI uses the same priors and likelihood as GIMs and a point estimate (BaCaDI*) and may be viewed as an ablation of GIMs that does not amortize the inference of M and I. For JCI-PC, we evaluate the classical (JCI-PC) (Brouillard et al., 2020) and feature vector variant (JCI-PC Context), which extends the graph with context nodes for the perturbation features. Since all causal approaches except GIMs are limited to predictions for seen perturbations 7, we also benchmark two block-box approaches: CondOT
Iiil functions (Gesztelyi
dom Hill functions, each assigned sampled target variable. The Hill functions define the shifts ab induced at their target variables I at given dosage e E R. This enables us to evaluate algorithms at two lev- els of generalization. In the partially distribiitii (OOD) task, we test on perturbations with new dosages
"Our code is available at https://gi com/ MoraSchneider/gin.
6

---

Generative Intervention Models for Causal Perturbation Modeling
A
C Intervention Effects
Linear Gaussian
Nonlinear Gaussian
Shift Intervention
SID (4)
Edge-F1 (t)
SID (1)
Edge-F1 ( (1)
(intervened ancestor)
1.0
1.0
* 0.0
200
0.5
200
0.5
0.5 1.0
0.0
0.0
No Intervention (intervened ancestor)
Partially
Fully 00D
Partially OOD
Fully 0OD
In Distribution
OOD
In Distribution
0.0 0.5
1.0
2.0
Target-F1(11
1.0
1.0
1.0
F1
1.0
1.0
1.0
Hard Intervention
(predicted intervened ancestor)
0.5
0.5
0.5
0.5
2.5
0.0
0.0
0.0
0.0
0.0
0.0
0.0 0.5 1.0
1.5 2.0
C
GIM
BaCaDi* GnlES
JCI-PC
Will JCI-PC (Context)
GIM
MLP
UT-IGSP
True
Data
Figure 2. Evaluating the learned causal structure and interventions. A & B: For linear and nonlincar SCMs, SID and Edge-F1 scores of inferred causal graphs (top) and Target-Fl scores (bottom) for in- and baselines at inferring the causal structure and intervention targets in nonlinear systems, while also generalizing to OOD perturbation features. C: Example intervention effects predicted by GIMs and MLPs compared to the ground truth (Hill function) on nonlinear SCMs. Data plotted at correspond the training data of 2 at dosage c. Lines show the predicted (mean) intervention effects at all dosages e [0, 2.5]. In top and center panels, GIMs accurately infer the true intervention effects, even though the training data suggests differently (true interventions target ancestors). In bottom panel, the GIM misses an intervention but instead predicts an intervention on ancestor.
perturbation features ry GIMs outperform all
(Bunne et al., 2022) and an MLP that predicts the mean perturbation shifts directly from y. Appendices D.4 and D.5 provide details on GIMs and the baselines, respectively.
results. Additional metrics and results for shift interventions and scale-free graphs are given in Appendix E.
Explicitly modeling the perturbation mechanisms im- proves structure learning nonlinear systems. Fig- ure 2 shows the accuracy of the estimated causal graphs compared the ground-truth SCM that generated the data. For linear systems (A), GIMs, BaCaDI", and both JCI-PC approaches predict similarly accurate graphs in terms of in- terventional implications (SID), while GnlES and UT-IGSP perform best. This aligns with previous findings that classi- cal methods may outperform continuous optimization-ban based approaches on linear, normalized data, suggesting the gap to GIMs may result from factors discussed by Reisach et al. (2021) and Ormaniec al. (2025). In nonlinear systems (B), GIMs significantly outperform the baselines across all met- rics, including BaCaDI". suggesting that the amortized infer- ence model can enhance causal discovery. JCI-PC (Context) performs on par with JCI-PC, indicating that merely in- cluding perturbation features without modeling their actual mechanisms is not sufficient for enhancing causal discovery.
5.3. Metrics
To evaluate the predicted perturbed distributions p(x;>" we need to compare various methods, some without an an- alytical p.d.f. To enable this, we directly compare samples X* from the predicted distribution with true samples X We report their catropy-regularizer Wasserstein distance (W2). the Euclidean distance between their means (Mean distance), and the average negative log- lig-liiel of the true data under a kernel density estimate of the predicted data (KDE-NLL). For the scRNA-seq data, we also report the Pearson correlation of the mean differential expression after perturbation as commonly used in scRNA-seg analysis. For the synthetic data, we also evaluate the inferred causal mech- anisms M and interventions Σ. Given estimated and true graphs G and G, we compute the structural intervention distance (SID) (Peters& & Bühlmann, 2015) and the Fl score of the edges (Edge-F1). Given predicted targets i, we report the Fl score given the true I (Target-FlI (see Appendix D.3).
GIMs generalize unseen perturbation features and ac- curately infer their mechanisms of action. Figure 2 also presents the Target-F1 scores for in-distribution, partially OOD, and fully OOD perturbations. GIMs significantly outperform all causal baselines in predicting intervention targets except for UT-IGSP, which performs best in linear settings. Moreover, the GIM framework is the only method
6. Results
In Figures 2 and 3, we present the results on synthetic data with hard interventions. Figure 4 analyzes the effect of the information content in y, and Figure reports the SciPlex3

---

Generative Intervention Models for Causal Perturbation Modeling
In Distribution
Partially OOD
Fully OOD
Mean Distance (
(4)
Mean Distance (4)
W2 (4)
Mean Distance (
30
I
20
4
20
2
10
10
GIM
BaCaDi
JCI-PC
UT-IGSP
MLP
CondOT
Observational
Figure 3. Benchmarking the predicted distribution shifts on nonlinear Gaussian SCMs. W2 and Mean Distance for SCMs under perturbations with: seen dosages and targets (left, in distribution), novel dosages with seen targets (center, partially OOD), and novel targets with dosages (right, fully OOD). Metrics medians all perturbations for given dataset. GIMs demonstrate robust predictive performance, matching the best causal baselines in distribution and performing on par with the best black-box approaches fully OOD Observational baseline uses the unperturbed distribution as prediction. Vertical line separates mechanistic and black-box approaches. Boxplots show medians and interquartile ranges (IQR). Whiskers extend to farthest point within 1.5 IQR from the boxes.
capable of making predictions for unseen perturbations. For partially OOD perturbations. GIMs continue identify true intervention targets with Target-F1 scores comparable to in-distribution settings. While lower accuracy in the fully OOD setting is expected, GIMs achieve test-time scores that match or even exceed the baselines' in-distribution perfor- mance.
tion on a variable in this scenario, it correctly predicts an intervention on the parent of that variable in the estimated causal graph.
GIMs predict distributions as accurately as black-box approaches fully OOD, In Figure 3, we report the ac- curacies of the predicted distributions of each approach, irrespective of whether the underlying model M matches the ground truth SCM. For in-distribution perturbations, the end-to-end predictions of GIM are competitive with the best causal baselines across all evaluated metrics. All methods appear to capture meaningful perturbation effects, as their performances surpass both the observational prediction and CondOT. Despite not explicitly modeling the causal mecha- nisms, the MLP consistently performs best at predicting the effects of seen perturbations across all metrics.
Figure C illustrates the perturbation mechanisms learned by GIMs and the MLP baselines for selected variables along- side the true Hill functions defining the interventions. These examples emphasize the differences between black-box pre- dictions and the GIM approach, which explicitly infers the mechanisms of action of a perturbation while allowing for end-to-end predictions. For instance, in the top and center panel, the samples of the four training perturbations suggest different variable shifts from those of the true underlying atomic intervention, because
of the ancestor variables is also intervened upon. While the MLP predicts shifts based on the observed data, the GIM correctly recovers the true Hill functions. In the bottom panel, GIM learns an incorrect intervention mechanism, yet it remains consistent with the estimated causal model. Although GIM misses an interven-
Only GIMs and black-box methods can make predictions on unseen perturbations. Partially OOD, GIMs generalize well and achieve scores comparable to those of the in-distribution setting. On fully OOD perturbations, all baselines perform worse. Yet, the predictions of GIMs and the MLPs are better than the observational data and thus nontrivial, improving significantly over CondOT. While the MLP performs best in partially OOD, GIMs tend to yield better results for fully OOD interventions, while providing a mechanistic explana- tion for their predictions via their internal causal model.
Varying Number of PCA Components (fully OOD)
0.6
%6 explained
30
Generalization requires that the perturbation features contain sufficient information. To allow for general- ization, the features 7 must allow for predictions about the (mechanistic) effects of unseen perturbations. To study this in GIMs, leverage our proposed PCA representation to encode the atomic interventions as features 7 (Section 5.1). Figures 4 and 11 show how the predictive performance and target identification of GIMs degrades as we vary the num- ber of PCA components contained in 7, keeping everything else fixed. When y contains sufficient information, GIMs effectively generalizes to completely unseen perturbations.
20
0.2
10
0.0
8
12
15
No PCA
unexplained variance of PCA
Observational
GIM
Figure 4. Predictive accuracy (W₂) of GIMs relative infor- mation contained in y. Fully OOD perturbations in nonlincar systems with hard atomic interventions. GIMs' predictive perfor- mance improves monotonically with the information content in
8

---

Generative Intervention Models for Causal Perturbation Modeling
Mean Distance (4)
Pearson Correlation (T)
CA12
TFPI2
18
1.0
175
Pomaii Taap
0.5
5
0.0
2
0
2
Gene Expression
Gene Expression
GIM MLP
CondOT
Naive
Observational
True data
Figure 5. Generalization to heldout drug-dosage combinations the SciPlex3 scRNA-seg drug perturbation data (Srivatsan et al., 2020). Mean Distance and Pearson Correlation between predicted and true data (left). Examples of true and predicted marginal distributions the 2 top-ranked marker genes for the Givinostat perturbation A549 cells (right). Naive corresponds predicting training data from the same drug and cell type at the closest lower dosage available in the training data (1,.M). GIMs achieve the most robust predictions, quantitatively and qualitatively, and the only mechanistic (causal) generative models applicable for unseen y
GIMs predict qualitatively accurate perturbation distri- butions in real-world settings. In Figure 5, we present the results on the scRNA-seq perturbation data. The predic- tions of GIMs and the MLP baseline are most accurate as measured by Mean Distance and correlation (see Figure 5, left; for KDE-NLL, see Figure 12 in Appendix E.4). Here, we do not report W2, GIMs predictions are continuous p.d.f.s, whereas all other baselines' predictions (including MLP) as well as the true data are quasi-discrete (transformed and normalized integer counts), rendering distributional dis- tances difficult to compare. Both GIMs and MLPs achieve better accuracy than naive approach that predicts the re- sponse of the next-lower dosage in the training data. When inspecting the marginals for selected genes (Figure 5, right), we find that GIMs capture the distribution shifts qualitatively most accurately, leveraging the full (zero-inflated) genera- tive model of the data. This shows that GIMs can provide meaningful predictions in settings with model mismatch.
Acknowledgements
This research was supported by the European Research Council (ERC) under the European Union's 's Horizon 2020 research and innovation program grant agreement no. 815943 and the Swiss National Science Foundation under NCCR Automation, grant agreement 5INT40 180545, We thank Scott Sussex for his valuable comments and feedback.
Impact Statement
Our work studies causal models, which are relevant to any field aiming at predicting the effect of interventions, pertur- bations, or policy changes. If algorithms are grounded in the causal structure of the true data-generating process, they may be less prone to unexpected behavior or wrong conclu- sions, provided the assumptions underlying our approach are verifiable in the given application. Since our work aims at ad advancing machine learning and statistics, there are many potential societal and scientific consequences of our work, none which we feel must be specifically highlighted here.
7. Conclusion
We introduced generative intervention models, general approach for using causal models to make predictions about perturbations of a data-generating process. GIMs learn to generate the atomic effects of a perturbation in a jointly- estimated causal model only based on its features. This allows us to tackle two key challenges in real-world decision- making: (1) inference of the causal mechanisms underlying Sata-generating process, which makes GIMs inherently interpretable, and (2) generalization beyond the observed perturbations, where existing causal inference methods are not applicable. This dual capability distinguishes GIMs from existing work and enables experts to incorporate prior knowledge about the model and intervention mechanisms. Future work may use our approach in fully Bayesian treat- ment (Lorch et al., 2021), to design optimally informative perturbations (Agrawal et al., 2019; Toth et al., 2022; Ailer et al., 2024), or scale it to larger systems (Lopez et al., 2022).
References
Agrawal, R., Squires, C., Yang, K., Shanmugam, K., and Uhler, C. ABCD-Strategy: Budgeted experimental de- sign for targeted causal structure discovery. In The 22nd International Conference on Artificial Intelligence and Statistics, pp. 3400-3409. PMLR, 2019.
Ailer, E., Dern, N., Hartford, J. S., and Kilbertus, N. Tar- geted sequential indirect experiment design. Advances in Neural Information Processing Systems, 38:122029- 122053, 2024.
Barabási, A.-L. and Albert, R. Emergence of scaling in random networks. Science, 286(5439) 1999.
Brouillard, P., Lachapelle, S., Lacoste, A., Lacoste-Julien S., and Drouin, A. Differentiable causal discovery from

---

Generative Intervention Models for Causal Perturbation Modeling
interventional data. Advances in Neural Information Pro- cessing Systems, 33:21865-21877, 2020.
Glorot, X. and Bengio, Y. Understanding the difficulty of training deep feedforward neural networks. In Pro- ceedings of the thirteenth International Conference on Artificial Intelligence and Statistics, pp. 249-256. JMLR Workshop and Conference Proceedings, 2010.
Bunne, C., Krause, A., and Cuturi, M. Supervised train- ing of conditional Monge maps. Advances in Neural Information Processing Systems, 35:6859-687 2022.
Gonzalez, G., Herath, I., Veselkov, K., Bronstein, M., and Zitnik, M. Combinatorial prediction of therapeutic pertur- bations using cousally-inspires neural networks. bioRxiv, 2024.
Bunne, C., Stark, S. G., Gut, G., Del Castillo, J. S., Levesque, M., Lehmann, K.-V., Pelkmans, L., Krause, A., and Rätsch, G. Learning single-cell perturbation sponses using neural optimal transport. Nature Methods, 20(11):179--1768, 2023.
Hägele, A., Rothfuss, J., Lorch, L.,, Somnath, V. R., Schölkopf, B., and Krause, A. BaCaDI: Bayesian causal discovery with unknown interventions. In International Conference on Artificial Intelligence and Statistics, pp. 1411-1436. PMLR, 2023.
Chickering, D. M. Optimal structure identification with greedy search. Journal of Machine Learning Research, 3 (Nov):507-554, 2002.
Chu, D., Zabet, N. R., and Mitavskiy, B. Models of transcrip- tion factor binding: sensitivity of activation functions to model assumptions. Journal of Theoretical Biology, 257
Hao, M., Gong, J., Zeng, X., Liu, C., Guo, Y., Cheng, X., Wang, T., Ma, J., Zhang, X., and Song, L. Large-scale foundation model on single-cell transcriptomics. Nature Methods, pp. 1-11, 2024.
(3):419-429, 2009.
Cui, H., Wang, C., Maan, H., Pang, K., Luo, F., Duan, N., and Wang, B. Segpt: toward building a foundation model for single-cell multi-omics using generative ai. Nature Methods, pp. 1-11, 2024.
Hauser, A. and Bühlmann, P. Characterization and greedy learning of interventional Markov equivalence classes of directed acyclic graphs. Journal of Machine Learning Research, 13479),2409-2464, 2012.
Cuturi, M. Sinkhorn distances: Lightspeed computation of optimal transport. Advances in Neural Information Processing Systems, 26, 2013.
Hetzel, L.,, Boehm, S., Kilbertus, N., Günnemann, S., Theis, F, et al. Predicting cellular responses to novel drug per- turbations at single-cell resolution. Advances in Neural Information Processing Systems, 35:26711-267221 2022.
Cuturi, M.,
L., Tian, Y., Bunne, C.,, Davis, G., and Teboul, O. Optimal transport tools (ott): A jax toolbox for all things wasserstein. arXiv preprint arXiv:2201. /2324, 2022.
Ibragimov, I. A. and Has' minskii, R. Z. Statistical Estima- tion Asymptotic Theory. Springer- Springer-Verlag. New York, 1981.
Dibaeinia, P. and Sinha, S. SERGIO: a single-cell expres- sion simulator guided by gene regulatory networks. Cell Systems, 11(3):252 22111,,,,111111 2020.
Jang, E., Gu, S., and Poole, Categorical reparameteriza- tion with gunbel-softmax. In International Conference Learning Representations, 2017.
Eberhardt, F., Glymour, C., and Scheines, R. On the number of experiments sufficient and in the worst case neces- sary to identify all causal relations among n variables. In Conference on Uncertainty in Artificial Intelligence, 2005.
Kamimoto, K., Stringa, B., Hoffmmann, C. M., Jindal, K., Solnica- -reieal, L., and Morris, S. A. Dissecting cell iden- tity via network inference and in silico gene perturbation. Nature, 6 614(7949):742-75,, 2023.
Erdōs, P.. Rényi, A., et al. On the evolution of random graphs. Publication of the Mathematical Institute of the Hungarian Academy of Sciences, 5(16:17-60, 1960.
Kingma, D. P. and Ba, J. Adam: A method for stochastic optimization. In International Conference on Learning Representations, 2015.
Gamella, J. L., Taeb, A., Heinze-Deml. C., and Bühlmann, P. Characterization and greedy learning of Gaussian struc- tural causal models under unknown interventions. arXiv preprint Xrxi:2221 14897, 2022.
Lee, H.-C., Danieletto, M., Miotto, R., Cherng, S. T., and Dudley, J. T. Scaling structural learning with NO-BEARS to infer causal transcriptome networks. In Pucific Sympo- sium on Biocomputing 2020, pp. 391 4022 World Scien- tific, 2019.
Gesztelyi, R., Zsuga, J., Kemeny- -ekke, A., Varga, B., Juhasz, B., and Tosaki, A. The Hill equation and the origin of quantitative pharmacology. Archive for History of Exact Sciences, 66:427-438, 2012.
Lopez, R., Regier, J., Cole, M. B., Jordan, M. L., and Yosef, N. Deep generative modeling for single-cell transcrip- tomics. Nature Methods, 15(12):1053-105 2018.
10

---

Generative Intervention Models for Causal Perturbation Modeling
Lopez, R., Hütter, J.-C., Pritchard, J., and Regev, A. Large- scale differentiable causal discovery of factor graphs. Ad- vances in Neural Information Processing Systems, 35: 19290-19303, 2022.
Peters, J., Janzing, D., and Schölkopf, B. Elements of causal inference: foundations and learning algorithms. The MIT Press, 2017.
Reisach, A., Seiler, C., and Weichwald, S. Beware of the simulated DAG! Causal discovery benchmarks may be easy game, Advances in Neural Information Processing
Lorch, L., Rothfuss, J., Schölkopf, B., and Krause, A. DiBS: Differentiable Bayesian structure learning. Advances in Neural Information Processing Systems, 34:24111- 24123, 2021.
Systems, 34:27772-27788 2021.
Roohani, Y., Huang. K., and Leskovec, J. Predicting tran- scriptional outcomes of novel multigene perturbations with GEARS. Nature Biotechnology, 42(6))927-935, 2024.
Lotfollahi, M., Wolf, F. A., and Theis, F. single-cell perturbation responses. Nature Methods, 16 (8):715-721, 2019.
segen predicts
Rothenhäusler, D., Heinze, C., Peters, J., and Meinshausen, N. BACKSHIFT: Learning causal cyclic graphs from unknown shift interventions. Advances in Neural Infor- mation Processing Systems, 28, 2015.
Lotfollahi, M., Klimovskaia Susmelj, A., De Donno, C.,, Hetzel, L., Ji, Y., Ibarra, L., Srivatsan, R., Naghipour- far, M., Daza, R. M., Martin, B., et al. Predicting cellular responses complex perturbations in high-throughput screens. Molecular Systems Biology, 19(6)e11517, 2023.
Sadybekov, A. V. and Katritch, V. Computational ap- proaches streamlining drug discovery. Nature, 616(7958); 673-685, 2023.
Maathuis, M. H., Colombo, D., Kalisch, M., and Bühlmann, P. Predicting causal effects in large-scale systems from observational data. Nature Methods, 7(4):247-244 248, 2010.
Schölkopf, B. Causality for Machine Learning, pp. 765-804. Association for Computing Machinery, New York. NY, USA, I edition, 2022. ISBN 9781450395861.
Maddison, C., Mnih, A., and Teh, Y. The concrete dis- tribution: A continuous relaxation of discrete random variables. In Proceedings of the International Conference on Learning Representations, 2017.
Schölkopf, B., Locatello, F., Bauer, S., Ke, N. R., Kalch- brenner, N., Goyal, A., and Bengio, Y. Toward causal representation learning. Proceedings of the IEEE, 109(5):
Mooij, J. M., Magliacane, S., and Claassen, T. Joint causal inference from multiple contexts. The Journal of Machine Learning Research, 2 2111):919-44026, 2020.
612-634, 2021.
Scott, D. W. Multivariate density estimation: theory, prac- tice, and visualization. John Wiley & Sons, 2015.
Ormaniec, W., Sussex, S., Lorch, L., Schölkopf, B., and Krause, A. Standardizing structural causal models. In The Thirteenth International Conference on Learning Representations, 2025.
Spirtes, P., Glymour, C., and Scheines, R. Causation, pre- diction, and search. MIT press, 2000.
Squires, C., Wang. Y., and Uhler, C. Permutation- based causal structure learning with unknown intervention tar- gets. In Peters, and Sontag, D. (eds.), Proceedings of the 36th Conference on Uncertainty in Artificial Intel- ligence (UAI), volume 124 of Proceedings of Machine Learning Research, pp. 1039-1048 1048. PMLR, 03-06 Aug 2020.
Parascandolo, G., Kilbertus, N., Rojas-Carulla, M., and Schölkopf. B. Learning independent causal mechanisms. In International Conference Machine Learning, pp. 4036 4046.4044. PMLR, 2018.
Pearl, J. Causality. Cambridge university press, 2009.
Srivatsan, S. R., McFaline-Fiiuurro. J. L., Ramani, V., Saun- ders, L., Cao, J., Packer, J., Pliner, H. A., Jackson, D. L., Daza, R. M., Christiansen, L., et al. Massively multi- plex chemical transcriptomics at single-cell resolution. Science, (67 2020.
Pervez, A., Locatello, F., and Gavves, S. Mechanistic neural networks for scientific machine learning. In Salakhut- dinov, R., Kolter, Z., Heller, K., Weller, A., Oliver, N., Scarlett, J., and Berkenkamp, F. (eds.), Proceedings of the 41st International Conference on Machine Learning, vol- ume 235 of Proceedings of Machine Learning Research, pp.
Tejada-Lapuerta, A., Bertin, P., Bauer, S., Alice, H., Bengio, Y., and Theis, J. Causal machine learning for single-cell genomics. Nature Genetics, pp. 1-12, 2025.
PMLR, 21-27 Jul 2024.
Peters, J. and Bühlmann, P. Structural intervention distance for evaluating causal graphs. Neural Computation, 27(3): 771-799, 2015.
Toth, C., Lorch, L., Knoll, C., Krause, A., Pernkopf, F,, Peharz, R., and Von Kügelgen, J. Active Bayesian causal

---

Generative Intervention Models for Causal Perturbation Modeling
inference. Advances in Neural Information Processing Systems, 35:16261-16275, 2022.
Vowels, M. J., Camgoz, N. C., and Bowden, R. D'ya like DAGs? A survey on structure learning and causal discov- ery. ACM Computing Surveys, 15(4):1-36, 2022.
Wang, Y., Solus, L., Yang, K., and Uhler, C. Permutation- based causal inference algorithms with interventions. Ad- vances in Neural Information Processing Systems, 30, 2017.
Yang. K., Katcoff, A., and Uhler, C. Characterizing and learning equivalence classes causal DAGs under in- terventions. In International Conference on Machine Learning, pp. 5541-5550. PMLR, 2018.
Yu, H. and Welch, J. D. Perturbnet predicts single-cell responses to unseen chemical and genetic perturbations. BioRxiv, pp. 2022-07, 2022.
Zheng, X., Aragam, B., Ravikumar, P. K., and Xing, E. P. DAGs with NO TEARS: Continuous optimization for structure learning. Advances in Neural Information Pro- cessing Systems, 31, 2018.
12

---

Generative Intervention Models for Causal Perturbation Modeling
A. Additional Background
Atomic interventions in a causal model refer to any external manipulations or changes that modify the variables or causal mechanisms within the network (Peters al., 2017: Vowels al., 2022). Formally. an intervention on a variable 2; corresponds to replacing the local conditional distribution pj in Equation (2) by a new interventional distribution pt, while the conditionals of the variables not subject the intervention remain unchanged (Schölkopf al., 2021; Schölkopf, 2022). Thus, interventions are sparse or local in nature, meaning that they only affect the targeted variables and do not directly influence the other variables of a system.
We typically distinguish between two types of atomic interventions: soft and hard interventions (Peters et al., 2017; Pearl, 2009; Eberhardt et al., 2005). Soft (or imperfect) interventions modify the conditional probability distributions of the intervened-upon variable without changing the structure the causal graph G. This means that the functional dependencies on the parents remain. One example for soft interventions are so-called shift interventions (Rothenhäusler al., 2015), which modify the conditional distributions of the target variables according to shift vector 2b (#11 ..., v/a) R. that
(1; = x + sb XG. 0, x/) = pr(e, = xG: 0).
(13)
Hard (or perfect) interventions, on the other hand, eliminate the dependencies of the intervened variables on their causal parents so that the interventional conditional simplifies to
P. (z : 0, = (z:
(14)
ultimately altering the structure of the causal graph. Real world examples from biology include perfect gene knockout experiments for hard interventions and gene silencing for soft intervention (Squires et al., 2020).
B. Derivations
B.1. Continuous Model Representation
Similar to Lorch et al. (2021), we derive an expression for the joint distribution of the model using the latent probabilistic graph model. We get
p(M", ΦD = p(Z0)p(0p(()p(())Iz)))) Z.000)
/ p(G DI 0, ØG
= =r(Z)p(0))(o)
=(())p(0)p(O) / |PG Z) II p p( (K) (. dG
k=1
/ p( Z)
K
/
= -(Z))(0)p())
II
| 10:00 p(() (1)) (k) |G,,, (1100 (M dG
k=1
B.2. Gradient Derivations
In the following, we detail the derivation of the gradients of the complete-dat log joint density, log loM &. D) (defined in Equation (9)), that allow us to use muddient-based optimization for the MAP estimation of the continuous representation of the causal model M., {Z, 0} and the parameters of the generative model Φ
Derivation of Vz Pog Φ
Vz logglo/M...ooo D)
Vologp(22 0) + Vz log 0)
= Vzlogp(Z) + Vzlog ()G])))p(01 G) dG
+ Vz log (fNGIZ)
K
II
((I(k) (1: (xx
G.O.I (h) dG
13

---

Generative Intervention Models for Causal Perturbation Modeling
Vz f p(G ()p(O) p () GG
G) dG Vz Z) Z) J p( (b:)(h))(())) G, G,0, dī dG
Vz (og
K fp(G Z) I ((Zik) (3;((())) ) (X( G, C, dz dG
During optimization, we sample i.dd G p G 2). ((m) p(0 Z) and [(m,k) p p( [ ()) for E {l,..,, MJ to upproximate the expectations. To ensure that we can compute the gradients through the samples, we use the reparame- terization trick when sampling 0 and vp, and the Gumbel-sigmoid trick when sampling G. This allows us to compute the gradients as follows:
Vz lolgoppMM, D)
1 Em Vzp(O(m) G(!))
1 EM
22 p(X G @(m), I)), Em=1 l (K G ((m), I(m,k])
Vzloog log + M
Σ ((0(m) ((m))
M
Derivation for Vo log Vologp(Ma. 0,
Vo loggp(Mar Φ D)
Vglogp(Maa) + (agloo(()MMar
K Vo log (FNGIIZ)O G) + Volog log (/NGIZZ z)
II / p(z(k)) 1 (t)) G, (G,0,I(0) d dG R=
Similar to computing the gradients with respect to Z, we sample from the respective distributions to approximate the expectations. use the reparameterization trick when sampling 0 and sb.
Derivation for V l ())
V (MM, Φ D) V log lologp(D) Ma. Φ)
V log + V log
(/
(
|
Z)
II k=1
K
/
p(I(h) (8) (xxx G, I()) az aG
During optimization, we sample from the respective distributions to approximate the expectations. We use the reparameteri- zation trick when sampling @ and &p and Gumbel-sigmoid trick when sampling I(A)
C. Computational Complexity
The computational complexity of GIMs is comparable to prior causal modeling approaches such as Hägele et al. (2023) and Brouillard et al. (2020) since we use similar causal model classes. The only additional overhead introduced by our new modeling framework is the forward pass through the GIM MLP, which does not affect the asymptotic complexity. Formally, computational complexity of learning GIMs is:
o (T [d (pz + HM* M Npower) + M Thotal LM H3 + K -LLGIM Ham + HOM (p+ + M .d))
where:
number of training steps,
d: number of variables in x,
pz: rank parameter of Z,
M: number of MC samples,
nowwer: number of power iterations used in NO-BEARS acyclicity,
n: total number of samples across environments
14

---

Generative Intervention Models for Causal Perturbation Modeling
LGIM, L depth of GIM MLPs and causal mechanisms, respectively,
HGIM, M width of GIM MLPs and causal mechanisms, respectively,
p: number of perturbation features in %
K: number of perturbation environments.
D. Details on Experimental Setup
In the following section, we provide details on our synthetic data generation, preprocessing the Sciplex dataset and implementation details on GIMs and all baseline approaches.
D.I. Synthetic Data Generation
linear and nonlinear Gaussian causal mechanisms as defined in Equation (7) and Equation (8), respectively. For linear causal mechanisms, linear functions model variable's dependency on its parents and thus the parameters O represent these regression coefficients. For data generation, we randomly sample the parameters & uniformly from [-3, -0.25] U [0.25, 3) to make sure they are bounded away from zero and exclude non-significant relationships. Further, we always use a constant variance for all variables 2;, thus exp(a;)² = exp(o)² = 0.1. Similarly, for nonlinear causal mechanisms, a variable's 's dependency on its causal parents defined via a nonlinear MLP. We thus have one MI.P for each variable, so we have in total d networks. We use MI.ss with one hidden layer with units and the tanh activation function. We sample the weights 0 from standard Gaussian with mean 0 and variance 1. Again, we use 2 = 0.1.
Causal Mechanisms For synthetic data generation,
We generate random atomic interventions, by first sampling the intervention targets I. Specifically, uniformly sample the targets from all possible subsets of variabels containing exactly one two elements. Next, we use randomly sampled Hill functions to compute the intervention parameters v/ acting
Atomic Interventions
a fixed target variable. For each mondomll-ssmppeed target variable I;, we separately sample A; ~ 4[-6, 22 U [2, 6] and define its Hill function as (c) = /1((1 + (0.75/c)*). Given the targets and their Hill functions, we then generate four distinct interventional environments for training by evaluating h,, (c) at the dosages c E {0.5, 1, 1.5, 2]. Akin to drug perturbations, each of these four environments shares the same targets but has different intervention effects sb, which are governed by the same Hill function evaluated at the different dosages c.
We consider types of atomic interventions: hard and shift interventions, which are detailed in Appendix A. For settings with hard interventions, the intervention parameters vb, which are specified via the Hill functions, describe the mean of the interventional conditional distribution for a targeted variable. Specifically, the interventional conditional of a targeted variable is defined as (z; (bb) = A (1:; vb exp(dintere) Contrary to the observational conditionals, we choose exp(Aintere) )² = 0.5 to ensure that the interventional distributions remain distinct from the observational ones even at small dosages, where the mean of the interventional conditional is close to zero, regardless whether a variable has parents the original graph. For settings with shift interventions, the intervention parameters v. which are specified via the Hill functions, merely shift the targeted variables without introducing additional noise as defined in Equation (13) (Rothenhäusler et al., 2015).
Synthetic Datasets The training dataset consists of total of 160 perturbational datasets (corresponding 40 Hill functions considered at 4 different dosages), where each one has n = 50 samples, and one observational dataset with no 800 samples. The partially -itto--distribttion test dataset consists of 200 perturbational contexts corresponding to the same 40 hill functions and intervention targets from the training dataset, but evaluated at 5 different dosages, e E (0.25, 0.75, 1.25, 1.75, 2.25]. Finally, the fully out-of-distribution test dataset consists of 80 perturbational contexts corresponding to 20 newly sampled hill functions, which are evaluated at dosages c E {0.5, 1,5,22)
D.2. SciPlex3 Data
Single-cell RNA sequencing data, like the SciPlex3 dataset (Srivatsan et al., 2020), usually comes with technical particulari- ties that require careful preprocessing. Thus, we apply standard preprocessing methods. First, we normalize the data such that each cell has the same total count, equal to the median count of the unperturbed data, We then filter cells to retain only
15

---

Generative Intervention Models for Causal Perturbation Modeling
those with at least 200 counts, and genes keep only those expressed in at least 20 cells. Next, we apply logarithmic transformation, ã log(1 x)- Since we focus subset of genes, we first identify the top 1000 highly variable genes, and subsequenly reduce this set to the top 50 marker genes. To determine these, we rank all genes and select those with the highest absolute scores.
We limit our analysis to four drugs (Belinostat, Dacinostat, Givinostat, and Quinostat) from the epigenetic regulation pathway. For each drug perturbation, we construct corresponding perturbation features by concatenating the one-hot encoding the drugs, one-hot encoding of the dosages, and rescaled dosage values, ensuring these values fall within [0. 1]. Specifically, we apply a hard-coded mapping of the dosages given by: 10n.M + 0.2, 100nAM + 0.4, 1μM 0.6, and 10p + 0.8M1 Thus, the perturbation vector 9-dimensional vector, with y {0, 1) x {0. x [0.
In total, we consider 12 different experimental settings, corresponding to three cell types and four drugs. Specifically, in each setting, we hold out one drug at its highest dosage (10μMM) as test dataset to evaluate generalization. The remaining perturbations are used for training, resulting in a training dataset x 3 + -3 + 1 16 perturbations.
D.3. Metries
Structural Interventional Distance (SID) We compare the predicted graph G to the true causal graph G using the SID (Peters & Bühlmann, 2015), which quantifies their closeness based on their implied interventional distributions. is defined as
SID(G, G) = # ( A(1.D)4
j E D
if E xG. Xa, is not valid adjustment set for (i,j) in G if j & xc.
}.
where D denotes the descendants of i in G. We refer to Peters & Bühlmann (2015) and Peters et al. (2017) for a detailed definition of valid adjustment sets and the SID.
Fi-score of predicted edges (Edge-F1). Treating the inference problem as classidication problem, we compute the Fl-score of the edges in predicted graph G with respect to the true causal graph G. The Fl-score is the harmonic mean of precision and recall, which is often prefered ove e.g. accuracy as the graphs are sparse. It is defined as
2TP
FI G) = 2TP +P+FN
where TP, FP and FN refer to the number of true positives, false positives and false negatives, respectively.
Fi-score of predicted intervention targets (Target-F1) We additionally report the Fi-score of the predicted intervention targets with respect to the true intervention targets (i, I), using the definition provided above.
Wasserstein distance W2 We compare the predicted samples X* to the true samples X using the W2 distance with a entropic regularization (Cuturi, 2013) provided in the OTT library (Cuturi et al., 2022) with
ΣΣ P X 8911 cHIP) 1/2 P (X.) 1=1
TT(X,XX") =
min
Here and m are the number of samples in X' and X', and U is the set of (n x m) transport matrices with U = (P e Raxm Pl. = 1/n. 1, and 1 = 1/m- n 1, with 1, denoting n-dimensional vector of ones. H(P) is the entropy with (PP) = E., Piy log P.j - 1. e is the regularization parameter and we always use c = 0.1.
Mean Distance We compare the predicted samples X to the true samples X using the Euclidean distance of their empirical means and is deñned
Mean-Distancee2 X" X*) = |1 1=
where ji" and 1²" are the empirical means of X* and X*, respectively.
16

---

Generative Intervention Models for Causal Perturbation Modeling
Average negative log-likelihood under kernel density estimation (KDE-NLL) We kernel density estimation (KDE) with Gaussian kernel and bandwidth chosen according to Scott's rule (Scott, 2015) obtain density based on the predictions X". Then, we evaluate the average negative log-likelihood of the true samples X based on the KDB-estimate
n
( KDE-NLL(X', X') == 72
log p(x;)
m
X
with p(x;) m bw
K
bw
where n and m are the number of samples in X* and X respectively, bw denotes the bandwidth and K the Gaussian Kernel.
Pearson Correlation between means We compare the predicted samples X to the true samples X using Pearson correlation on their means. Specifically, the Pearson correlation is defined
Pearson Correlation X* X*) Ei (it N:) 15²)
VE (A) P0S 1T
Here we denote ji" and μ* as the empirical means of X* and X respectively, and f and 2rt are the average values of ji and μ², respectively.
D.4. GIM Implementation Details
In the following, we detail the implementation of GIMs used for our experiments reported in Section 6.
Causal Mechanisms M: On synthetic data, we choose the causal inference model to have no mismatch with the generated data. Specifically, for linear data, the linear Gaussian inference model defined in Equation (7), and for nonlinear data, we use the nonlinear Gaussian inference model defined in Equation (8), with the MLPs containing single hidden layer of 5 nodes and tanh-activation.
Because single-cell gene expression data is subject to dropout effects and remains positive even after our transformation, we model it with a zero-inflated log-normal (log.V) inference model, which satisfies log(z) ~ W. The corresponding conditional for each variable v; 2 0 is given by
p(a; xG. ; 0) nT 8(x,) +
(o,loooNW,, (G, x)), exp((.))
(15)
where 5(x) is Dirac delta at zero, and n, the learned probability of zero-inflation and o, represents the learned noise scale for each variable r, Following the approach for nonlinear Gaussian mechanisms, we model the causal mechanisms with an MLP, selecting its architecture through hyperparameter tuning and using tanh activation function.
In all experiments, we set the variance in the model priors in Equation (12) to n = - (where is the number of system variables), n = 0.1, n 4. In the prior of the GIM, we choose n = 0.1.
In our proposed approach, we ensure the acyclicity G via the model prior p(M.) (see Equation (12)). Specifically, we use the continuous acyclicity constrained NO-BEARS proposed by Lee et al. (2019), who show that the spectral radius, i.c. the largest absolute eigenvalue, of adjacency matrix is 0 if and only if the corresponding graph is acyelic. As suggested by Lee et al. (2019), we use power iteration compute and differentiate through the largest absolute eigenvalue:
u G
cGG) = P(G)
TT T
where we update for t = 30 steps according to the following update rule:
u - 2 T G//llu" TG|2u + Gu/[Gulle
and u, e E IR are initialized randomly.
17

---

Generative Intervention Models for Causal Perturbation Modeling
GIM
We use the function Pe map the perturbation features to the probabilities of the Bernoulli distributions over the intervention targets I. For this function, we consistently apply MLP with two hidden layers, each containing 100 nodes, and tanh activation function.
GIMs allows us to model different intervention types. On synthetic data, we model interventions in each setting according to the true underlying intervention type. For shift interventions, 2/ represents the shift parameters, and h maps the perturbation features to the means of a Gaussian distribution over vp. For hard interventions, &b includes both the interventional means and the interventional noise scale for a targeted variable. Here, hg consists of two separate MLPs: one mapping to the interventional means and the other the interventional noise scale, both using the same architecture. For all synthetic settings, the MLPs for the intervention parameters & have 2 hidden layers with 100 nodes and tanh activation function. On drug perturbation data, we assume hard atomic interventions. Thus, v/ represents all parameters of the zero-inflated log- Gaussian conditional, that is the zero-inflation probability, the means and the noise scale. Thus he consists of three separate MIPs, one for each parameter. We choose the architecture of the MLPs based on hyperparameter tuning and use a tanh-activation function.
Optimization Details and Initialization We employ gradient- basedd optimization to obtain MAP estimates for Ma. and d. In Appendix B.2 we provided the gradients of the posterior allowing us to use Adam optimization (Kingma & Ba, 2015) with a learning rate of 0.001. On synthetic data, we use 30000 steps and on drug perturbation data, we use 100000 steps. For the Monte Carlo approximations, we use a sample size of TMc = 128. We apply a cosine annealing schedule to the coefficient, 3z. which controls the sparsity of the intervention targets.
We initialize Z by sampling from a standard Gaussian distribution. initial parameters via Omit (O, 0.1 I). For nonlinear mechanisms, use the Glorot normal (Glorot & Bengio, 2010) to initialize the parameters mitt The log- variances are initialized with zero, 0, 0. Last, use the Glorot initialization for the weights of GIM, ø. This initialization procedure is important to avoid inducing a bias at the start of the optimization process.
When using linear causal mechanisms, we sample the
We initialize the Lagrangian multiplier and the penalty coefficients with A = and po = 1 respectively. Every 100 optimization steps, we evaluate based on 100 holdout samples if the objective converged. If it converges and the acyclicity constraint, cGG) not satisfied, we update the Lagrangian multiplier and penalty according to
A411 = d + pe c(G))
(16)
1411 = 2/4--
(17)
D.5. Baselines
BaCaDi (Hägele et al., 2023): BaCaDi is causal discovery approach for settings with unknown interventions and uses pradeent-based variational inference. Thus it lears posterior over the causal model and the intervention targets and parameters of each perturbation context. In order to obtain a point estimate, we implement only one particle in BaCaDi's inference procedure. Additionally, we adjust the inference model and its priors to match GIM's inference model. To distinguish our implementation from the original one, we call our implementation BaCaDi* Through the modifications, B acts as a counterpart to GIM, which uses a similar inference model but does not model the perturbation mechanisms and thus cannot generalize to unseen perturbations.
GnIES (Gamella et al., 2022): GnIES, an extension of the GES algorithm (Chickering, 2002), is score-based causal discovery approach for settings with unknown interventions. GnIES optimizes penalized likelihood score by alternately optimizing the set of intervention targets and the graph equivalence class. It assumes that the noise distribution of targeted variables changes across interventional contexts. As it assumes a linear model, there a model mismatch for nonlinear systems. GnlES only estimates a graph equivalence class. Thus, for computing Edge-F1, we favor them by correctly orienting undirected edges when present in the ground truth graph. In case of a wrongly predicted edge, we orient randomly. Since GnlES outputs set of intervention targets for the entire dataset instead of perturbation-wise intervention targets, end-to-end perturbation predictions are not straightforward.
JCI-PC (Mooij et al., 2020) The Joint Causal Inference (JCI) framework (Mooij et al., 2020) extends general causal discovery methods to accommodate data from various interventional contexts, including those with unknown interventions.
18

---

Generative Intervention Models for Causal Perturbation Modeling
The authors suggest including information describing system's context as additional nodes in the causal graph. Subse- quently, standard causal discovery algorithms can be applied to the augmented graph. We implement this approach using the PC algorithm (Spirtes et al., 2000), which relies sequence of conditional independence tests to infer equivalence class graph estimate. We use a Gaussian conditional independence test for all settings, since the KCI-test did not terminate in reasonable time. We implement two types of JCI-PC: First, we add one node for each perturbational context (similar to the experiments conducted in (Hägele et al., 2023; Brouillard et al., 2020)). Second, we add the perturbation features as additional nodes to the causal graph (JCI-PC (Context)) Importantly, this approach does not infer the intervention targets for each perturbation context. As JCI-PC only estimates a graph equivalence class, we compute causal discovery metries similar to GnlES. In order evaluate the end-to-end perturbation predictions, sample graph from the equivalence class and use linear Gaussian SCM with maximum likelihood parameter and variance estimates as the learned model (Hauser & Bühlmann, 2012).
UT-IGSP (Squires et al., 2020): UT-IGSP extends the GSP (Yang et al., 2018; Wang et al., 2017) algorithm to settings with unknown intervention targets. Specifically, UT-IGSP learns a graph equivalence class and intervention targets by optimizing regularized local scores based on local conditional independence tests and permutation search. We compute causal discovery metrics similar to our GnIES implementation, as UTI-GSP also only estimates a graph equivalence class. To evaluate the end-to-end predictive performance, we follow the same procedure as for JCI-PC.
CondOT (Bunne et al., 2022): CondOT learns a global optimal transport map conditioned on perturbation's 's context variable, that effectively maps measure onto another. This global map generalizes and ultimately allows predictions for new contexts by conditioning them. For our settings, we use the perturbation features y as the context variable and map the distribution of observational data D(0) to the different perturbational ones D(R), As the data of our experiment is already low-dimensional, we apply CondOT in the original dataspace.
make
Multilayer Perceptron (MLP): We use an MLP to predict the shift of the perturbational distribution (A) relative to the observational one X. from Yk- Specifically, for each perturbational context, we calculate the shifts = pa 44o where /4 and (i are the means of the datasets (A) and X respectively. The MLP is trained predict these shifts from the perturbation features using a mean-squared error loss. When predicting the effects of new perturbations we shift the observational data X by the predicted shift, resulting in X {a; + 8" Ma
D.6. Hyperparameter Tuning
In our experiments on synthetic data, we benchmark GIM in different generative settings (linear/ nonlinear, hard/ shift interventions). For each setting, we perform separate search of hyperparameters using 10 additionally generated datasets and choose the hyperparameters corresponding to the optional average KDE-NLL on the test datasets. Since the causal baselines are limited to predictions for seen interventions, we add the test dataset to the training data. Table 1 shows the range of hyperparameters investigated.
For our experiments on the Sciplex3 data, we use similar approach. From each training dataset, we randomly hold out one drug at the extreme dosage and then select the hyperparameters to optimize the average KDE-NLL of the held-out samples. The search ranges are provided in Table 2
E. Additional Experimental Results
F.1. Additional results on Erdosrenyi Graphs
In Figure 6, we report an additional metric for measuring the accuracy of the predicted distribution, NLL-KDE, on nonlinear Gaussian SCMs with ER graphs, so the same settings reported already in Figure 3. Additionally, we show in Figure 6 the predictive accuracy for linear Gaussian SCMs.
In Figure 7, we report additional results for Gaussian SCMs with ER graphs, but with shift, atomic interventions underlying the perturbations. These results mostly match the results for hard interventions. We select hyperparameters described above.
19

---

Generative Intervention Models for Causal Perturbation Modeling
Table I. Hyperparameter tuning for the experiments on synthetic data Section 6.
Method
Hyperparameter
Search Range
GIM
Sparsity graph
8M & (50, 100, 500}
Sparsity intervention targets
8r E {1, 10, 50, 100)
Inverse temperature sigmoid
TE (0.01, 1,
JCI-PC
Significance level
JCI-PC (Context)
Significance level
GnlES
Penalization parameter
E [BIC penalizaationn,,5,,,111111 |
UT-IGSP
Significance level (CI)
acr E 6, le 5, 11--6-1666,,61-44000000000000000
Significance level (invariance)
E {le-6, le-5, 1e- 000---_000___00--000000000000000
BaCaDi
Sparsity graph
8M E {50, 500
Sparsity intervention targets
87 {1, 10, 50,50,0000000 100}
Inverse temperature sigmoid
T {0.01, 0.1, 10)
MLP
Hidden layers
E ([50), (50, 50J, (100), [100, 100), (100, 100, 100))
Learning rate
E {0.001, 0.01}
Number of iterations
Rater E {1000, 5000, 20000, 40000)
CondOT
Hidden layers
E [64], [64, 64), [128, 128), [64, 64, 64, 64]}
Learning rate
Ir 6 [0.001, 0.01}
Number of iterations
Rater {60000, 80000, 100000}
Table 2. Hyperparameter tuning for the experiments on SciPlex3 data in Section
Method
Hyperparameter
Search Range
GIM
Inference Model: Hidden layers
Riny C {[20), [20, 20), [20, 20, 20])
GIM: Hidden layers
RGIM E {1100, [100, 100), [100, 100, 100]}
Sparsity graph
BM E 1,
Sparsity intervention targets
87 E {0.1, 10}
MLP
Hidden layers
E ([50), (50, 50), [100), [100, 100). [100, 100, 100])
Learning rate
I {0.001, 0.01)
Number of iterations
{1000, (1000,5000, 10000, 20000, 40000)
CondOT
Hidden layers
E [[44). (64, 64), [128, 128), [64, 64, 64, 64])
Learning rate
Ir E {0.001, 0.01}
Number of iterations
Rater {60000, 80000, 100000}
E.2. Scale-Free Graphs
In Figure 8 and Figure we report experimental results for synthetic data with scale-free graphs using hard atomic and shift atomic interventions, respectively. The experimental setup the same as described in Section 5 with the only difference being that use the hyperparameters that were optimized for settings with ER graphs for GIM and all other baselines. The results are mostly consistent with the ones on ER graphs and align with findings discussed in Section 6.
E.3. Impact of Hyperparameter Selection and Feature Information Model Performance
As discussed in Appendix D.6, select the hyperparameters to optimize the predictive accuracy of the perturbed distribution, as measured by KDE-NLL. The accuracy of the estimated causal structure M is expected to be closely related to the predictive accuracy. Yet, full identification is not strictly necessary for good predictions of the perturbed density under y. To put in perspective, we report in Figure 10 the causal discovery metrics for GIM and all causal baselines, when selecting the hyperparameters based on Edge-F1. GIM and all causal basclines achieve higher accuracy compared to the setting where hyperparameters are selected based on KDE-NIL. However, in many real-world settings, the true causal model is unknown
20

---

Generative Intervention Models for Causal Perturbation Modeling
In Distr KDE-NIL 149
Partially O0D
Fully GOD KDE-NLL (4)
KDI-MLL (4)
25
25
20
15
MLP
CondOT
Oeservational
(a) Nonlinear Mechanisms.
In Distribution
Partially OOD
Fully OOD
W;
Mean Distance
KDE KDE. (1)
W; (4)
Mean Distance |KDE-NLL
W,
Mean Distance (4 JLL 44)
4
25
15
30
20
10
15
5
15
GIM
JCI-PC
UT-GSSP
BaCaDi
MLP
CondOT
Observational
(b) Linear Mechanisms
Figure 6. Benchmarking the predicted distribution shifts Gaussian SCMs with ER graphs and hard interventions. (a) KDF-NLI for SCMs with nonlinear mechanisms. (b) W2. Mean Distance and KDE-NLL NLL for SCMs with linear mechanisms. We show perturbations with training dosages and targets (left, in distribution), novel dosages with targets (center, partially OOD), and novel targets with seen dosages (right, fully OOD). Metrics are medians over all perturbations for a given dataset. Observational baseline uses the unperturbed distribution prediction. Vertical line separates mechanistic and black-box approaches. Boxplots show medians and interquartile ranges (IQR). Whiskers extend to the farthest data point within 1.5 IQR from the boxes.
and thus the Edge-F1 is not accessible for selecting the hyperparameters. Further, the goal of our work is to use causal model to ultimately predict the perturbed distribution and thus choosing hyperparameters based on metrics measuring the predictive accuracy might be more suited.
In Figure 11 we additionally report the accuracy of the inferred intervention targets when varying the number of principal components used for constructing the perturbation features. While for in-distribution perturbations, the performance of GIM converges with about principal components, more information in is required for generalizing.
E.4. Additional results for SciPlex3 (Srivatsan al., 2020)
In Figure 12, we report the KDE-NLL for the predicted distribution of GIM and bascline approaches on SciPlex 3 data. Additionally, in Figure 13, show the marginal histograms of the GIM and MLP prediction in one selected perturbation scenario. This visual example shows how the MLP merely shifts the observational distribution, while GIM, using the zero-inflated model, qualitatively better matches the true perturbations distributions.
21

---

Generative Intervention Models for Causal Perturbation Modeling
Linear Gaussian (ER)
Nonlinear Gaussian (ER)
SID 4)
Edge-f1 11
SID (4)
Edge-F1 (1)
300
300
0.75
0.50
100
100
0.25
0.00
GIM
BaCaDi*
KI-PC
JCI-PC (Context)
UT-IGSP
GnlES
(a) Causal Discovery
Linear Gaussian (ER)
Nonlinear Gaussian (ER)
Distribution
Portially 0oD
Fully DOD
In Distribution
Partially 00D
Rully 000
1.0
1.0
1.0
1.0
1.0
1.0
Target:11 [1
0.0
0.0
0.0
GIM
BoCaDh
KCPPC
UT-IGSP
(b) Target Identification
In Distribution
Partially OOD
Fully OOD
W;
(
(1)
(4)
W.(4)
30
15
30
3
20
20
20
2
15
JCI-PC
UT-IGSP
BaCaDie
CondOTT
Observational
(c) Predictive Distribution for Linear Mechanisms.
In Distribution Distance
Partially OOD
Fully OOD
(1)
[1]
Distance JKDE-NLL
Wy [
Mean
(L LIDE-NLL (1)
30
30
25
30
20
25
20
1.0
2
20
20
25
10
35
20
GIM
JCH-PC
UTHIGPP
BaCaDi
MLP
CondOT
Observational
(d) Predictive Distribution for Nonlinear Mechanisms.
Figure 7. Evaluating the learned causal structure, interventions and distribution shifts Gaussian SCMs with ER graphs and shift interventions. (a) SID and Edge-Fl scores of inferred causal graphs for linear (left) and nonlincar (right) Gaussian SCMs. (b) Target-Fl scores in- (d) Wz. Mean Distance and KDE-NLL for SCMs with linear (c) and nonlinear (d) mechanisms. We show perturbations with training dosages and targets (left, distribution), novel dosages with seen targets (center, partially OOD), and novel targets with seen dosages (right, fully OOD). Metrics are medians over all perturbations for given dataset. Observational baseline uses the unperturbed distribution as a prediction. Vertical line separates mechanistic and black-bo box approaches. Boxplots show medians and interquartile ranges (IQR). Whiskers extend the farthest data point within 1.5 IQR from the boxes.
disttiibititnionnn perturbation features y for linear(left) and nonlinear (right) Gaussian SCMs. (c) and
22

---

Generative Intervention Models for Causal Perturbation Modeling
Linear Gaussian (SF)
Nonlinear Gaussian (SF)
SID 1)
Edge-F1
SID
())
Edge-F1 (1)
0.75
300
200
200
0.50
GIM
BaCaDi
(ClnneCC)
GnlES
(a) Causal Discovery
Gaussian (SF)
Nonlinear
(SF)
In Distribution
Partially 0oD
Fully DOD
In Distribution
Partially OOD
Fully 000
1.0
1.0
1.0
1.0
1.0
1.0
Target:Fl
0.0
GIM
BoCaD
JC-PC
UT-IG5P
(h) Target Identification
In Distribution
Partially OOD
Fully OOD
(4)
(4)
(1)
ILDE-NLL
W, (4)
(4
15
30
25
20
1.0
20
15
25
20
10
20
15
JCI-PC
Observational
(c) Predictive Distribution for Linear Mechanisms.
In Distribution
Partially OOD
Fully OOD
W3 (
(
W5 (4)
30
25
20
20
15
20
JCI-PC
CondOTT
Observational
(d) Predictive Distribution for Nonlinear Mechanisms.
Figure 8. Evaluating the learned causal structure, interventions and distribution shifts Gaussian SCMs with SF graphs and hard interventions. (a) SID and Edge-F1 scores of inferred causal graphs for linear (left) and nonlinear (right) Gaussian SCMs. (b) Target-FI scores for in- and -ut-o--distribution perturbation features y for linear(left) and nonlinear (right) Gaussian SCMs. (c) and (d) W2. Mean Distance and KDE-NLL for SCMs with linear (c) and nonlinear (d) mechanisms. We show perturbations with training dosages targets (left, distribution), novel dosages with seen targets (center, partially OOD), and novel targets with seen dosages (right, fully OOD). Metrics are medians over all perturbations for a given dataset. Observational baseline uses the unperturbed distribution as prediction. Vertical line separates mechanistic and black-box approaches. Boxplots show medians and interquartile ranges (IQR). Whiskers extend the farthest data point within 1.5 IQR from the boxes.
23

---

Generative Intervention Models for Causal Perturbation Modeling
Linear Gaussian (SF)
Nonlinear Gaussian (SF)
SID
Edge-F1 [1)
EdgeF1
0.75
0.75
200
0.50
200
0.50
0.25
0.25
100
0
0.00
0.00
GIM
BaCaDi
JCH-C
Wele JCI-CC (Context)
UT-IGSP
GnIES
(a) Causal Discovery
Linear Gaussian (SF)
Nonlinear Gaussian [SF)
Distribution
Partially
Fully
Partially
Fully
Target-F1
0.5
0.5
0.5
0.0
0.0
0.0
0.0
0.0
0.0
GIM
BaCaDi*
JCHPC
UT-I6SP
(b) Target Identification
In Distribution Distance
Partially OOD
Fully OOD
Mean Distance JKDE-NLL
Maan (4 NLL
30
30
25
2
15
20
10
20
15
Observational
(c) Predictive Distribution for linear mechanisms.
In Distribution Distance
Partially OOD
Fully OOD
KDE
Mean Distance |XDE--NLL
Mean
30
30
20
25
3
25
20
30
2
20
3
(d) Predictive Distribution for nonlincar mechanisms.
Figure 9. Evaluating the learned causal structure, interventions and distribution shifts on Gaussian SCMs with SF graphs and shift interventions. (a) SID and Edge-FI scores inferred causal graphs for linear (left) and nonlinear (right) Gaussian SCMs. (b) Target-Fl scores for in- and out-of-distribunt perturbation features y for lincar(left) and nonlinear (right) Gaussian SCMs. (c) and (d) W2, Mean Distance and KDE-NLL for SCMs with linear (c) and nonlinear (d) mechanisms. dosages and targets (left, distribution), novel dosages with targets (center, partially OOD), and novel targets with dosages (right, fully OOD). Metrics are medians over all perturbations for a given dataset. Observational baseline uses the unperturbed distribution prediction. Vertical line separates mechanistic and black-box approaches. Boxplots show medians and interquartile ranges (IQR). Whiskers extend to the farthest data point within 1.5 IQR, from the boxes.
show perturbations with training
(ER)
Nonlinear Gaussian
SID (1
Edge-F1 11
SID
Edge-F1 (11
300
1.0
300
200
200
0.50
100
100
0.25
GIM
BaCaDi*
JCI-PC
222
JCI-PC (Conte(()n)eet)
UT-ISSP
GnIES
Figure Causal Discovery when Selecting Hyperparameters based on Edge-F1. SID Edge-Fl scores inferred causal graphs for linear (left) and nonlinear (right) Gaussian SCMs with hard, atomic interventions. Boxplots show medians and interquartile ranges (IQR). Whiskers extend to the farthest data point within 1.5 TQR from boxes.
24

---

Generative Intervention Models for Causal Perturbation Modeling
Varying Number of PCA Components
Distribution
Fully Out of Distribution
Target-11 (1)
0.8
1.0
explained
0.8
1.0
%% % exsplinna
0.8
0.8
Target:F1
0.6
0.6
0.6
0.6
0.4
0.2
BaCaDi*
GIM
BaCaDi
explained variance of PCA
Figure Target Identification (Target-F1) GIMs relative information contained in y. Target-Fl for in-distribution (Ieft) and fully OOD perturbations (right) in nonlinear systems with hard atomic interventions. Boxplots show medians and interquartile ranges (IQR), Whiskers extend to farthest data point within 1.5 IQR from the boxes.
KDE-NLL (4)
125
120
70
60
50
GIM MLP
CondOT
Observational
Naive
Figure 12. KDE-NLL for heldout drug-dosage combinations on SciPlex3 scRNA-seq drug perturbation data (Srivatsan et al., 2020).
CA22
THH2
HNFIAAS5
IGPL2A481
ALDHIAL
KATB2
COL25A2
RBFOX3
GAPCH
AKhi8i
AC080197.1
CBACAM22I
ACSL3
GHKSA
MUSSAC
KHTLS
MK67
CHOALNACTI
AGR2
UNC002506
NPLLL
CENF
UNC00552
DHR52
SUC2A3
ACTO1
HNOCSA
LNC522004
NORFALL
TOP2A
True
GIM
MLP
Observational data
Figure 13. True and predicted marginal distributions on SciPlex3 data. 50 top-ranked marker genes for the Givinostat perturbation on A549-cells.
25


Reviews:
[{"review_id": "ZpGLGRHcXS", "reviewer": "Reviewer_jgn6", "summary": "This paper studies the problem of causal perturbation modeling to recover the causal structure and intervention targets given perturbation features from several interventional environments. The use-case in this paper is the gene perturbations in the biology domain. The authors propose generative intervention model (GIM) which learns to map perturbation features to atomic intervention vectors (1-sparse) in a jointly-estimated causal model. This modeling enables the ability to infer the effect of unseen perturbation features during inference. Thus, this approach is robust to distribution shifts in the feature perturbation space. Experiments are conducted on synthetic and scRNA-seq drug perturbation data to show the effectiveness of the proposed method.\n\n## Update After Rebuttal\nI thank the authors for the detailed response to my questions and concerns. I believe this is an interesting work that relaxes the major assumption of having intervention targets available in causal representation learning. Furthermore, I believe the application in biology is quite interesting. Therefore, I lean towards **acceptance** of this paper.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed comments and feedback. We are glad to hear that you find our approach addresses an important problem in causal generative modeling with an interesting use-case with great potential. You raised several important points that we address in detail below. \n\n> overall objective can use more clarity. \n\nWe revised the manuscript in lines 95ff, 133ff, 154ff to further clarify.\n\n> What is the main difference between this approach for extrapolating to unseen perturbations and the causal representation learning method proposed by Zhang et al?\n\nThank you for pointing us to the work by Zhang et al. (2023). Their setting differs from ours: their causal model is defined over latent (unobserved) variables, whereas GIM models interventions over observed variables (e.g., gene expression). This makes GIM’s predictions more interpretable, as it directly models causal effects between semantically meaningful variables. While it is unknown in both approaches how the perturbation modifies the underlying causal model, the two methods tackle this challenge in fundamentally different ways. Zhang et al. (2023) infer interventions for each perturbation in the training data, independent of any features. In contrast, we propose to learn how interventions are induced by the observable perturbation features $\\gamma$. Since Zhang et al. (2023) do not incorporate such features, their model *cannot generalize to novel perturbations* beyond those seen during training, similar to Hägele et al. (2023). Zhang et al.’s approach is limited to combinations of previously observed interventions, while GIMs enable predictions for entirely unseen new perturbations by conditioning explicitly on $\\gamma$.\nWe also note that, because GIMs model the causal variables directly, known identifiability results can be extended to GIMs (see our response to Reviewer 8P9S). Overall, we view the two approaches as complementary, exploring different but important aspects of causal modeling with perturbations.\n> meaning of the perturbation vector\n\n\n$\\gamma$ encodes any observable information about the perturbation applied in a given environment, such as drug properties and dosage in a drug perturbation experiment. For example, in the scRNA-seq setting, we define it as a one-hot encoding of the drug along with its dosage as the perturbation features. The perturbation vectors $\\gamma$ indeed differ across experimental environments. For each environment (i.e., context or experiment), we have one perturbation vector $\\gamma \\in \\mathbb{R}^p$, and a corresponding collection of samples of the observed variables $\\mathbf{x}$, such as gene expression levels. In line 78, we adjusted the manuscript to clarify the meaning of $\\gamma$ and that it differs across environments. \n\n> this setting does not explicitly require access to interventional data. Rather, one requires only a perturbation vector. Is this interpretation correct?\n\nYes, your interpretation is correct. GIMs do not require access to interventional data in the sense of having knowledge of the atomic intervention $\\mathcal{I}$ corresponding to a perturbation — that is, which variables in the causal model were directly manipulated and how. During training, we assume access to $K$ different perturbation experiments, where for each experimental context we observe perturbation features $\\gamma$, along with samples of observed variables $\\mathbf{x}$ collected under that condition. At test time, GIMs only require a perturbation vector $\\gamma^*$ to predict the system’s response. This setup reflects many realistic biological settings, where intervention targets are unknown but perturbation metadata is available.\n\n> multiple interventions can also be supported by simply inferring each one separately from the perturbation features. However, how does this work in practice?\n\nOur approach handles multiple perturbations (i.e., environments) by inferring a distribution over interventions $\\mathcal{I}$ — that is, intervention targets and parameters — for each perturbation individually based on its perturbation features $\\gamma$. In practice, the generative intervention model consists of two neural networks, $g_\\phi$​ and $h_\\phi$​, which map $\\gamma$ to the parameters of this distribution: $g_\\phi$​ outputs Bernoulli probabilities over targets, and $h_\\phi$​ outputs the mean and variance of a Gaussian over intervention parameters. Crucially, the parameters $\\phi$ are shared across all perturbation environments. Given a different $\\gamma^{(k)}$, the same model predicts a distinct intervention distribution for that specific perturbation. Combined with the causal model $\\mathcal{M}$, this yields a full predictive distribution over the system variables $\\mathbf{x}$.\nIn short, multiple perturbations are handled by applying the shared generative model to each perturbation feature vector $\\gamma$, enabling per-environment inference in a unified way.", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "Yes", "theoretical_claims": "N/A", "experimental_designs_or_analyses": "Yes, I checked the metrics used, including the Structural Intervention Distance (SID) and Wasserstein distance for evaluation. Furthermore, I checked the soundness of the data-generation with respect to the structural causal models and the empirical setup.", "relation_to_broader_scientific_literature": "Overall, this paper tackles an important problem in causal generative modeling. Specifically, for domains such as biology, predicting the effects of unseen interventions in critical for scientific discovery. Along a line of work that focuses on causal models for biological applications, such as intervention design and causal representation learning, this work explores a different approach to interventional inference.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "## Strengths\n\n- Overall, the paper is written well with clear intuitions about the problem of causal intervention modeling. Furthermore, this paper seems to make weaker assumptions about access to interventional data, specifically that one would only have access to the perturbation features and not the entirety of the interventional data for all environments, which is a more realistic scenario in practice.\n- The application in single-cell gene perturbations is certainly an interesting use-case with great potential impact, especially for reasoning about hypothetical perturbations outside the support of the training distribution.\n- The partial OOD (novel dosages with seen targets) and fully OOD (novel targets with seen dosages) is a robust way to evaluate the effects of unseen perturbations. Specifically, the fully OOD setting shows promising results compared to other methods.\n\n## Weaknesses\n\n- The overall objective can use more clarity. The idea seems to be to parameterize a network $\\phi$ to infer the intervention target given the perturbation vector. To learn a joint causal model, one objective is to learn a latent variable $z$ to model the causal graph $G$ probabilistically. The overall goal is to maximize a log-likelihood with respect to the causal model and intervention prediction network over the dataset.", "other_comments_or_suggestions": "N/A", "questions_for_authors": "- What is the main difference between this approach for extrapolating to unseen perturbations and the causal representation learning method proposed by Zhang et al? In this method, the authors show identifiability guarantees from soft interventions and recover the causal factors and their relationships up to an equivalence class given sufficient interventional data.\n- Intuitively, what is the meaning of the perturbation vector? Is this simply the perturbation that generates K different interventional environments? Is it a value vector that takes a different value across environments?\n- From my understanding, this setting does not explicitly require access to interventional data. Rather, one requires only a perturbation vector. Is this interpretation correct?\n- The authors claim that multiple interventions can also be supported by simply inferring each one separately from the perturbation features. However, how does this work in practice?\n\nZhang et al. Identifiability Guarantees for Causal Disentanglement from Soft Interventions. NeurIPS 2023.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Ngy3xFBnZH", "reviewer": "Reviewer_8P9S", "summary": "The authors considered the problem of predicting the impact of interventions with applications in gene perturbation prediction. In particular, in some applications, when an intervention is performed, it is unknown which variables are intervened on. However, some features of the intervention might be known. The authors motivated this by an example in genomics, where features of a drug might be known but it is not clear what is the causal effect of a drug on regulatory pathways in a cell. The authors considered a causal modeling approach so that the outputs are easily interpretable. More specifically, they trained a generative model to map the features of a perturbation to a distribution over atomic interventions in the system. The experimental results showed that for the task of causal discovery or intervention target identification, the proposed method has better performance compared to previous work. Moreover, for predicting the effect of perturbation in out-of-distribution settings, the proposed method achieved similar performance compared to methods that do not use causal modeling in their architectures.\n\n### Update after rebuttal\nI increased my score to 3 after the discussion on the theoretical guarantees. However, as I mentioned before, the current paper is somehow a natural extension of Lorch et al. (2021) and Hägele et al. (2023), particularly in scenarios where intervention targets are not given or where features are not used to locate intervention targets. Therefore, I am not sure about the technical novelty in the current work and not giving a higher score.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed feedback. We are glad to hear that you find our problem setting relevant, the experiments comprehensive, and the interpretability of our approach valuable. You raised several important points that we address below. Please let us know if you have any remaining concerns or questions, and we would be happy to clarify further.\n> no theoretical guarantee [...] I should emphasize that this comment is also applied to some of the previous work.\n\nWe propose a new modeling framework that generalizes to unseen perturbations and demonstrates strong empirical performance on simulated and real-world data. While a full theoretical analysis is beyond the scope of this paper, we agree that identification guarantees are important. Below, we describe how existing results apply directly to GIMs, which allow contextualizing our results. We added these explanations to line 180. \nUnder standard assumptions, the MAP optimizers in GIMs identify the true intervention targets of the training perturbations, $\\mathbb{I}$, and a graph $\\mathbb{I}$-Markov equivalent to the true one. To establish this, we extend Theorem 2 from Brouillard et al. (2022). They prove that, under standard assumptions and with sufficiently small regularization, the score $\\mathcal{S}(\\mathbf{G},\\mathbb{I})$ is maximized by the true intervention targets and an $\\mathbb{I}$-equivalent graph. Assuming the GIM prior distributions have global support, we recover the same maximizers in the large-sample limit, because the posterior is dominated by the likelihood: $\\arg\\max_{\\mathbf{G}, \\mathbb{I}} \\mathcal{S}(\\mathbf{G},\\mathbb{I})=\\arg\\max_{\\mathbf{G},\\mathbb{I}}\\sup_{θ}\\log p(θ, \\mathbf{G},\\mathbb{I}\\mid\\mathcal{D})$. Assuming that there exists a mapping from features to the true targets and it can be expressed by $\\phi$, we have: $\\max_{\\mathbf{G}, \\mathbb{I}}\\sup_{θ}\\log p(θ,\\mathbf{G},\\mathbb{I}\\mid\\mathcal{D})=\\max_{\\mathbf{G},\\phi} \\sup_{θ}\\log p(θ,\\mathbf{G},\\phi\\mid\\mathcal{D})$. The MAP optimizers recover $\\mathbb{I}$ and a graph $\\mathbb{I}$-Markov equivalent to the true one. Identifiability of interventions for unseen perturbations depends on the informativeness of $γ$, and we empirically show the impact on predictive performance in Figure 4.\n> design choices/approximations are considered [...] why it is fine to consider them.\n[....] explain more about the design choice for target sparsity\n\nWe approximate expectations in the prior and likelihood using Monte Carlo sampling, and we use the Gumbel-Softmax relaxation for differentiable sampling of discrete variables. The MC estimators are consistent, while the Gumbel-Softmax introduces a bias that vanishes as the sigmoid temperature $τ \\to 0$. In practice, we fix $τ$. These approximations are widely used in the structure learning literature.\nOur prior design follows standard causal modeling assumptions, including sparsity and acyclicity for the graph and Gaussian priors over parameters for regularization. For the intervention targets, we use an L1-based sparsity-inducing prior, reflecting the assumption that typically only a few variables are affected per perturbation. This aligns with the sparse mechanism shift hypothesis, which suggests that interventions tend to induce localized changes (Schölkopf, 2022). All priors have full support; thus, in the large-sample limit, their impact diminishes.\n> what is $p$ [...]? What are $z_{0i}$ and $z_{1j}$ [...]?\n\n$p$ is a parameter controlling the rank of the matrices of $p(\\mathbf{G}\\mid Z)$. For $p\\geq d$, this distribution can represent any adjacency matrix without self-loops. In all experiments, we set $p=d$. The tensor $Z$ consists of two $d \\times p$ matrices; $z_{0i}$ and $z_{1j}$ refer to the $i$-th row of the first and the $j$-th row of the second matrix, respectively. We added this explanation in line 193.\n> [In Figure 2C, the] proposed method misses an intervention on a variable but predicts an intervention on its parent\n\nWe show this example to discuss how GIMs behave in practice and how they may infer consistent but nevertheless incorrect interventions in certain cases. Here, the true intervention targets $X_i$, and in the true graph, $X_j$ is a child of $X_i$. The intervention shifts the distribution of $X_i$, which affects $X_j$ via $p(X_j \\mid X_i)$. In the estimated graph, the edge is reversed and the model instead attributes the shift to an intervention on $X_j$. While this is incorrect under the true graph, it is consistent with the predicted graph and still describes the observed distribution.\n> main contributions compared to previous work\n\nPlease see our response to reviewer feSa.\n> method highly depends on $M$\n\nThe predictive performance can be affected by the number of MC samples $M$.  In all experiments, we used $M=128$. We additionally tested smaller $M$ in the nonlinear setting (Figure 3) and found no significant effect on accuracy:\n| |M=2|M=32|M=128|\n|---|---|---|---|\n|$W_2$|17.78 ± 5.35|17.69 ± 4.99|18.13 ± 4.33|", "claims_and_evidence": "The work is more experimental although it also provides proof for the derivations of some of the equations in the appendix. It seems that the authors made a good comparison with previous work.", "methods_and_evaluation_criteria": "They considered several metrics for comparison such as SID and Edge-F1 for causal discovery or $W_2$ and mean distance for the task of predicting the impact of perturbations. Moreover, they considered both synthetic and real datasets in their evaluations.", "theoretical_claims": "The work is experimental and there is no theoretical claim in the paper.", "experimental_designs_or_analyses": "I did not check the codes but did read the experiment section. Based on the text, it seems that the experimental results are sound and the authors also provided some explanations for the plots.", "relation_to_broader_scientific_literature": "One of the main applications of causal inference/discovery is in biology such as learning causal structures in gene regulatory networks or predicting the results of gene perturbations. This work is aligned with this line of research and aims to predict the effect of intervention when it is unknown which variables are intervened by giving the treatment.", "essential_references_not_discussed": "As far as I checked the related work, the authors cited the most relevant previous work.", "other_strengths_and_weaknesses": "Strengths:\n- The authors proposed a method to train a causal generative model that is more interpretable (for instance, checking which variables are intervened) compared to previous work, especially the unstructured models. \n\n- They performed comprehensive experiments in various settings, showing that the proposed method is better or on par with previous works.\n\nWeaknesses:\n- The main weakness is that there is no theoretical guarantee in the work. That being said, I should emphasize that this comment is also applied to some of the previous work. \n\n- In deriving the objective function of the optimization problem, several design choices/approximations are considered. It is not clear what are the impacts of such approximations and why it is fine to consider them.", "other_comments_or_suggestions": "It would be nice to study whether some theoretical guarantees can be added to the work in some specific settings such as linear settings. Moreover, it is good to characterize how the contained information in perturbation features $\\gamma$ affects perturbation prediction. For this purpose, some simple settings such as bivariate models might be considered to facilitate the analysis.", "questions_for_authors": "- In line 184, what is $p$ in the definition of the domain of $Z$? What are $z_{0i}$ and $z_{1j}$ in line 190? Please explain in more detail the modeling in line 190?\n- In eq. (11), please explain more about the design choice for target sparsity.\n- In the experiments, it was observed that the proposed method misses an intervention on a variable but predicts an intervention on its parent or ancestor. It is good to explain this observed phenomenon.\n- It seems that parts of the proposed method were borrowed from previous work such as Lorch et al. (2021) or Hagele et al. (2023). It is good to mention that other than the generative process of $I$ given $\\gamma$, what are the main contributions compared to previous work?\n- I think the performance of the method highly depends on $M$. It is good to study experimentally the effect of $M$ on the performance. I am geussing $M$ should be very large in order to have reasonable performance.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "TzELj5ZXaK", "reviewer": "Reviewer_feSa", "summary": "This paper studies the problem of predicting the unseen perturbation effect in a causal model. The authors propose the Generative Intervention Model (GIM) framework to learn the relationship between perturbation and distribution shift in a causal model. It is claimed that the GIM can predict perturbation features that do not show up in the training data. Detailed experiments on synthetic data and scRNA-seq are done to showcase the good performance of the GIM methods compared to other causal inference methods.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed comments and feedback. We are glad to hear that you find our work to be supported by clear and convincing evidence, as well as our experiments to be well designed. You raised several important points that we address in detail below.\n> Jointly estimating both the causal model and the generative intervention mapping seems to be a complex optimization problem, which may be computationally expansive.\n\nThe computational complexity indeed increases compared to settings, where only the causal model and intervention targets are directly inferred (e.g. BaCaDi (Hägele et al., 2023)). However, in practice, we find the overhead to be manageable. For example, in the 20-node nonlinear setting (as used for Figure 2B) and training on a GPU, BaCaDi* takes on average 11.56 ± 0.11 minutes, while GIM takes 20.71 ± 0.33 minutes. Both results are averaged over 20 random datasets. \n> add a concrete real-world example [....] & explain the formulation more clearly? In particular, what are known and what are unknown?\n\nA concrete real-world example is drug perturbation experiments, such as those on scRNA-seq data that we present in Section 6, Figure 5. We perform $K$ experiments, where in each one a drug is applied to cells and we measure their responses - typically via its gene expression profiles $\\mathbf{x}$. In experiment $k$, we observe the response of $n_k$ cells, yielding $n_k$ samples of $\\mathbf{x}$, denoted as $\\mathbf{X}^{(k)}$. The perturbation features $\\gamma^{(k)}$ capture observable drug properties, such as identity, molecular characteristics, or dosage. The full dataset $\\mathcal{D}$ consists of $K$ pairs of drug features and responses: $(\\gamma^{(k)}, \\mathbf{X}^{(k)})$ for $k = 1, \\dots, K$. Our goal is to predict the response of a new drug $\\gamma^*$, i.e. to predict the gene expression under this new drug.\n\nThe *latent* causal model $\\mathcal{M}$ defines the data-generating process of $\\mathbf{x}$. A perturbation $\\gamma^{(k)}$ influences this system through an *unobserved* intervention $\\mathcal{I}^{(k)}$. The distribution over these interventions, conditioned on $\\gamma^{(k)}$, is parameterized by *latent* parameters $\\phi$, which are shared across all experiments and enable generalization to new perturbations.\n\nIn summary, we observe $K$ pairs of $(\\gamma^{(k)}, \\mathbf{X}^{(k)})$. The causal model $\\mathcal{M}$, the interventions $\\mathcal{I}^{(k)}$, and the parameters $\\phi$ are unobserved and jointly inferred. We added this clarification in line 134. \n> similarity and difference between your work and [Lorch et al. (2021), Hägele et al. (2023)].\n\nGIMs introduce a new modeling framework that uses causal models to predict the effects of previously unseen, general perturbations, a setting that causal discovery methods like Lorch et al. (2021) and Hägele et al. (2023), are not designed to handle. While we adopt a similar prior design for the causal model $\\mathcal{M}$ and use gradient-based inference, we introduce a fundamentally different approach to modeling perturbations. \nLorch et al. assume that the intervention $\\mathcal{I}$ is observed for each perturbation; unlike GIMs, their method cannot be applied when interventions are unknown. Hägele et al. infer interventions for each training perturbation, but treat each perturbation independently, without leveraging any features. As a result, their model cannot generalize to novel perturbations beyond those seen during training. In contrast, we propose to learn how interventions are induced by the observable perturbation features $\\gamma$. Specifically, GIMs learn a shared generative mechanism that maps observable perturbation features $\\gamma$ to a distribution over latent interventions, enabling generalization to entirely new, unseen perturbations. Thus, while prior approaches focus on recovering the causal graph or inferring training interventions, GIMs are designed for the task of predicting the effects of new perturbations.\n\n> how acyclicity is ensured when sampling from the prior?\n\nWe enforce acyclicity via constrained optimization using the augmented Lagrangian approach, as introduced in Zheng et al. (2018) and also used in related works, such as Brouillard et al. (2020). While this encourages acyclic graphs during training, there is no hard guarantee that the estimated graph is cycle-free. In practice, we randomly break cycles if they occur.\n\n> provide the code in a later version\n\nWe will provide the full code in the camera-ready version.", "claims_and_evidence": "Claims made in the submission supported by clear and convincing evidence.", "methods_and_evaluation_criteria": "The paper uses synthetic data and scRNA-seq datasets in the experiments. The entropy-regularized Wasserstein distance between predicted and ground-truth distribution and the Euclidean distance between the mean are used as metric.", "theoretical_claims": "I did not check the derivations in Appendix B carefully, but they seem correct.", "experimental_designs_or_analyses": "The experiments are well-designed. Implementation details are provided.", "relation_to_broader_scientific_literature": "Several related works are mentioned in the Introduction part", "essential_references_not_discussed": "Relevant works are discussed.", "other_strengths_and_weaknesses": "Strengths: The proposed GIMs leverage perturbation features to predict distribution shifts in unseen conditions, which enables to model to generalize beyond seen distribution. Experimental results show that GIMs outperform classical approaches in identifying causal structures and intervention targets in nonlinear systems.\n\nWeaknesses: Jointly estimating both the causal model and the generative intervention mapping seems to be a complex optimization problem, which may be computationally expansive.", "other_comments_or_suggestions": "It would be good if the authors could add a concrete real-world example in the third section to illustrate the terminologies. It took me a while to understand all the notations. In particular, explanations of $\\gamma$ and how it may influence intervention are needed. Besides, i am still a bit confused about what is in the dataset $\\mathcal{D}$. Is intervention $I$ in the dataset?", "questions_for_authors": "Thanks for the great work! It will be great if the author can clarify the following questions. \n\n1. As I mentioned above, can you explain the formulation more clearly? In particular, what are known and what are unknown?\n\n2. The authors also mention in the paper that [1] uses a similar Bayesian method for structure learning. Can you explain the similarity and difference between your work and [1]? \n\n3. Can the author explain a bit about how acyclicity is ensured when sampling from the prior? \n\nI also wonder if the authors can provide the code in a later version. \n\n[1] Lorch, Lars, et al. \"Dibs: Differentiable bayesian structure learning.\" Advances in Neural Information Processing Systems 34 (2021): 24111-24123.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Nz75WxdVG9", "reviewer": "Reviewer_rdKL", "summary": "This paper presents Generative Intervention Models (GIMs), a causal modeling framework designed to predict the effects of perturbations in complex systems with unknown underlying mechanisms. GIMs establish a mapping between perturbation features and a distribution over atomic interventions within a jointly inferred causal model. This method facilitates the prediction of distribution shifts caused by unseen perturbations. The authors assess GIMs using synthetic datasets and single-cell RNA sequencing (scRNA-seq) drug perturbation data, demonstrating that GIMs outperform traditional causal inference methods in structure and intervention target identification while maintaining predictive accuracy comparable to black-box machine learning models.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed comments and feedback. We appreciate that you recognize the strong empirical results of the work and the potential impact in scientific domains such as biology.  You raised several important points that we address in detail below. Please let us know if you have any remaining concerns or questions, and we would be happy to clarify further.\n\n> it relies on black-box modeling and lacks theoretical guarantees or bounds on the correctness of its predictions. \n\nWhile our approach leverages neural networks, it is specifically designed *not* to be a purely black-box predictor. GIMs are built on a causal modeling framework with interpretable, mechanistic components: a structural causal model $\\mathcal{M}$ and interventions $\\mathcal{I}$, which model how perturbations manifest in $\\mathcal{M}$. This allows us to not only predict the effect of a perturbation, but also understand qualitatively how it alters the system’s underlying causal structure and mechanisms.\n\nWhile a full theoretical analysis is beyond the current scope, existing identifiability results from causal discovery with unknown interventions (e.g., Brouillard et al., 2020) apply to our setting under standard assumptions. We discuss this in detail in our response to Reviewer 8P9S. \n\n> technical novelty is limited, as it adopts an existing Bayesian structure learning framework, e.g., Brouillard et al.\n\nGIMs introduce a new modeling framework that uses causal models to predict the effects of previously unseen, general perturbations, a task that standard causal discovery methods are not designed to handle. While we build on existing techniques for inferring the causal model $\\mathcal{M}$ (e.g., prior design and gradient-based optimization), we introduce an explicit  model of how perturbation features induce interventions in the causal model. Specifically, GIMs learn a shared generative mechanism that maps observable perturbation features $\\gamma$ to distributions over latent interventions, enabling generalization to entirely new perturbations. In contrast, methods such as Brouillard et al. (2020) identify the unknown intervention targets for each training perturbation without incorporating their features, typically with the goal of recovering the causal graph. This limits their use to scenarios where the goal is recovering the causal graph and predicting effects for known perturbations, making them unsuitable for predicting the effects of new, unseen perturbations.\n\n>  What is the computational complexity of the proposed method?\n\nThe computational complexity of GIMs is comparable to prior causal modeling approaches such as Hägele et al. (2023) and Brouillard et al. (2020) since we use similar causal model classes. The only additional overhead introduced by our framework is the forward pass through the GIM MLP, which does not affect the asymptotic complexity. In terms of wall-clock time, for a nonlinear system with 20 variables (as in Figure 2B), training a full GIM takes on average 20.71 ± 0.33 minutes on a GPU, averaged over 20 random datasets. For reference, BaCaDi* takes 11.56 ± 0.11 minutes under the same conditions.\n\nFormally, the computational complexity of learning GIMs is:\n$$ \\mathcal{O}\\left( T\\cdot \\left[ d^2 \\cdot \\left(p_z + H_{\\mathcal{M}} + n_{\\text{MC}} n_{\\text{power}}\\right) + n_{\\text{MC}} \\cdot n_{\\text{total}} \\cdot d \\cdot L_{\\mathcal{M}} \\cdot H_{\\mathcal{M}}^2 + K \\cdot \\left( L_{\\text{GIM}} \\cdot H_{\\text{GIM}}^2 + H_{\\text{GIM}} \\cdot (p + d) + n_{\\text{MC}} \\cdot d \\right) \\right]\n\\right), $$ \nwhere:\n- $T$: number of training steps,\n- $d$: number of variables in $\\mathbf{x}$,\n- $p_z$: rank-parameter of $\\mathbf{Z}$,\n- $n_{MC}$: number of MC samples,\n- $n_{\\text{power}}$: number of power iterations used in NO-BEARS acyclicity,\n- $n_{\\text{total}}$: total number of samples across environments,\n- $L_{GIM}, L_{\\mathcal{M}}$: depth of GIM MLPs and causal mechanisms, respectively,\n- $H_{GIM}, H_{\\mathcal{M}}$: width of GIM MLPs and causal mechanisms, respectively,\n- $p$: dimension of $\\gamma$,\n- $K$: number of perturbation environments.\n\nWe added this to the supplementary material of our manuscript. \n\n> method applicable to complete gene datasets containing thousands of genes?\n\nLike comparable causal modeling approaches, GIMs do not scale easily to very large systems. Since our focus is on introducing a new modeling framework for perturbation effects, scalability is not a primary concern – in particular, because our experiments on scRNA-seq data demonstrate that GIMs can achieve strong predictive performance in real-world settings, even when applied to a subset of genes. That said, we agree that scalability is important and may be improved in future work. Approaches developed for scaling causal models, for example modeling the graph $\\mathbf{G}$ as a low-rank factor graph (Lopez et al., 2022), can be naturally integrated into our framework by modifying the causal model $\\mathcal{M}$ accordingly.", "claims_and_evidence": "The effectiveness of the proposed method is empirically well supported. \n\nHowever, it relies on black-box modeling and lacks theoretical guarantees or bounds on the correctness of its predictions. This limitation may restrict its applicability in real-world settings.", "methods_and_evaluation_criteria": "The authors conducted evaluations across various settings, including synthetic and scRNA-seq datasets. The evaluation criteria are well-chosen and provide detailed insights into the method's effectiveness.  \n\nOne experimental limitation is the lack of analysis on larger-scale settings with hundreds of thousands of nodes. (Authors subsample 50 genes for scRNA-seq experiments.) It raises the question of whether the method faces computational challenges when applied at such a scale.", "theoretical_claims": "This paper does not contain any theoretical claims for their method.", "experimental_designs_or_analyses": "The experimental design is sound. Some limitations regarding experiments are described in the evaluation section above.", "relation_to_broader_scientific_literature": "The key contributions of this paper are applicable to a broad range of scientific domains, particularly in the biological field.", "essential_references_not_discussed": "Overall, the paper presents key literature on causal modeling. However, since I am not familiar with generative intervention modeling, I am unable to assess whether the paper includes adequate baselines for comparison.", "other_strengths_and_weaknesses": "Strengths:  \n- Provides a comprehensive overview of existing literature. \n\nWeaknesses: \n- The technical novelty is limited, as it adopts an existing Bayesian structure learning framework, e.g., Brouillard et al., \"Differentiable Causal Discovery from Interventional Data,\" NeurIPS 2020.", "other_comments_or_suggestions": ".", "questions_for_authors": "- What is the computational complexity of the proposed method?  \n- Is the method applicable to complete gene datasets containing thousands of genes?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper investigates causal perturbation modeling to recover the causal structure and intervention targets, given perturbation features from several interventional environments. In particular, the authors proposed a generative intervention model (GIM) that learns to map these perturbation features to distributions over atomic interventions in a jointly-estimated causal mode. In addition, the authors applied the proposed method to gene perturbation prediction and obtained promising results. All the four reviewers agreed that the paper is well written, and the effectiveness of the proposed method is well supported by empirical results. Reviewers raised some concerns about assumptions, technical details, computational complexity, novelty, etc. Many of them have been addressed by the detailed responses from authors, but some reviewers still concern about the novelty of this work, compared to existing work such as Lorch et al. (2021) and Hägele et al. (2023).
Decision: Accept (poster)