Interpretable Causal Representation Learning for Biological Data in the Pathway Space
Authors: Jesus de la Fuente Cedeño
Robert Lehmann
Carlos Ruiz-Arenas
Jan Voges
Irene Marín-Goñi
Xabier Martinez de Morentin
David Gomez-Cabrero
Idoia Ochoa
Jesper Tegnér
Vincenzo Lagani
Mikel Hernaez
Venue: ICLR 2025 Poster (2025)
Keywords: Causal Representation Learning
Intepretability
VAE
Genomic Perturbations
Health
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)

Abstract: Predicting the impact of genomic and drug perturbations in cellular function is crucial for understanding gene functions and drug effects, ultimately leading to improved therapies. To this end, Causal Representation Learning (CRL) constitutes one of the most promising approaches, as it aims to identify the latent factors that causally govern biological systems, thus facilitating the prediction of the effect of unseen perturbations. Yet, current CRL methods fail in reconciling their principled latent representations with known biological processes, leading to models that are not interpretable. To address this major issue, in this work we present SENA-discrepancy-VAE, a model based on the recently proposed CRL method discrepancy-VAE, that produces representations where each latent factor can be interpreted as the (linear) combination of the activity of a (learned) set of biological processes. To this extent, we present an encoder, SENA-$\delta$, that efficiently compute and map biological processes' activity levels to the latent causal factors. We show that SENA-discrepancy-VAE achieves predictive performances on unseen combinations of interventions that are comparable with its original, non-interpretable counterpart, while inferring causal latent factors that are biologically meaningful.

Published as a conference paper at ICLR 2025
INTERPRETABLE CAUSAL REPRESENTATION LEARNING FOR BIOLOGICAL DATA IN THE PATHWAY SPACE
Jesus de la Fuente Robert Lehmann', Carlos Ruiz-Arenas' Jan Voges' Irene Marin-Goni> Xabier Martinez-de-Morentinn David Cobez-Cabreroos Idoia Ochoa Jesper Tegner", Vincenzo Lagani and Mikel Hernaez 6,+
CIMA University of Navarra, CCUN, IdiSNA, Pamplona, Spain.
TECNUN, University of Navarra, San Sebastián, Spain.
3 Biological and Environmental Science and Engineering Division, King Abdullah
University of Science and Technology (KAUST), Thuwal, Saudi Arabia
4 Dept. of Molecular Pharmacology and Experimental Therapeutics, Mayo Clinic, MN, USA
5 Institute of Chemical Biology, Ilia State University. Tbilisi 0162, Georgia
6 Center for Data Science (DATAI), University of Navarra, 31008, Pamplona, Spain.
Corresponding authors.
ABSTRACT
Predicting the impact of genomic and drug perturbations cellular function crucial for understanding gene functions and drug effects, ultimately leading improved therapies. To this end, Causal Representation Learning (CRL) constitutes of the most promising approaches, as it aims to identify the latent factors that causally govern biological systems, thus facilitating the prediction of the effect unseen perturbations. Yet, current CRL methods fail in reconciling their principled latent representations with known biological processes, leading to models that are not interpretable. To address this major issue, present SENA-discrepancy- VAE, a model based on the recently proposed CRL method discrepancy- VAE, that produces representations where each latent factor can be interpreted as the (linear) combination of the activity of a (learned) set of biological processes. To this extent, we present an encoder, SENA-6, that efficiently compute and map biological processes' activity levels to the latent causal factors. We show that -ENN--iicrepancc- VAE achicves predictive performances on unseen combinations of interventions that are comparable with its original, -on-interpretable counterpart, while inferring causal latent factors that are biologically meaningful.
INTRODUCTION
Causal Representation Learning (CRL) has raised in recent time promising approach for identify- ing the latent factors that causally govern the systems under study (Schölkopf et al., 2021; Ahuja et al., 2023). Among other disciplines, CRL have been recently applied on biological systems, providing precise testable predictions on causal factors associated with disease or treatment resistance (Zhang al., 2024; Lopez al., 2023). These methods usually operate on mixture of observational and interventional biological data, exploiting the distributional shift caused by the interventions with the goal of retrieving the causal latent factors and, possibly, how they mutually interact with each other. Concomitantly, Perturb-seq (Dixit et al., 2016) data have emerged as an ideal testbed for these type of analyses. This technology allows the gene expression profiling of single cells both in their unperturbed state and when
genes made functionally inoperative (c.g., through CRISPR knock-outs (KO) (Gilbert et al., 2014)). While generating expression profiles for thousands of cells across variety of experimental conditions is indeed advantageous in the CRL context, the high dimensionality of Perturb-seq data presents notable challenges for these models.
Deep learning approaches have allowed to predict transcriptional outcomes of novel (combinations of) perturbations in Perturb-seq data (Roohani et al., 2024; Cui et al., 2024; Gaudelet et al., 2024), or of known perturbations novel cell types (Lotfollahi et al., 2019). However, the best of our

---

Published as a conference paper at ICLR 2025
knowledge, there are only two works so far applying CRL to Perturb-seg data (Lopez al., 2023; Zhang et al., 2024), and we argue that these models severely lack interpretability, as the reconstructed latent factors cannot be directly reconciled with known biological processes, yielding talent factors that are difficult to interpret. Attempts made so far to boost interpretability in these models include computing associations between the reconstructed latent factors and the activity of known processes (Lopez et al., 2023), or arbitrarily selecting genes as representative of each latent causal factor (Zhang al., 2024). As also suggested in recent review (Tejada-Lapuert) et al., 2023), using biological processes as prior knowledge during the reconstruction process, rather than afterwards, may indeed boost the interpretability of the resulting models.
Related work. In the context of interpretable Representational Learning (RL), recent years have seen extensive applications of variational autoencoders (VAEs) to single-cell applications (Lopez et al., 2018). Some of these approaches enable interpretable latent factors either by enforcing gene-cell correspondence during training (Choi et al., 2023), by performing pathway enrichment analysis on linear gene embeddings (Zhao et al., 2021), or by modifying the VAE architecture to mirror user-provided gene-pathway maps (Seninge et al., 2021; Gut et al., 2021; Lotfollahi et al., 2023; Niyakan et al., 2024; Ruiz-Arenas et al., 2024a). Importantly, while these latter methods apply architectural changes that are similar in spirit to the ones we propose in this work, none of them present strong theoretical guarantees for the causal interpretation of their embeddings.
Therefore, in this paper we show how CRL algorithms can be extended in order to employ bio- logical processes (BPs) as prior knowledge, improving the interpretability of the resulting latent factors. We base our results on the CRL framework first introduced by Ahuja et al. (2023) and then expanded by Zhang et al. (2024). In particular, we present LEENA-discrepancy based on the recently proposed discrepancy-VAE. We show that SENA-discrepance VAE yields predictive performances comparable to the ones of its original counterpart on unseen combination of perturbations, while providing a mapping between latent factors and biological processes. To this end, we modify the discrepancy- VAF's encoder architecture (Zhang et al., 2024) and embed it with biological processes as prior knowledge. To our knowledge, this is the first effort to reconcile CRL with biological interpretability, achieving both principled identifiability and interpretability of causal latent factors in the biological pathway space.
CRL model
PRELIMINARIES AND BACKGROUND
2.1 CASUAL REPRESENTATION LEARNING
In what follows, we use the notation of Zhang et al. (2024), where we further use upper-case to denote random variables, lower-case to denote (inferred/observed) realizations of the random variables, upper-case bold to denote matrices, and lower-case bold to denote vectors. Let's assume that samples E R" are generated according to process governed by a set of latent variables U E Rd where d <<<>> n. These latent variables are not required to be independent from each other. Instead, each latent factor U, may be regulated by a subset of other latent factors, namely its parents Pa(U.), according to a structural mechanism U, + s; (;(PUU.)) Z.), where Z; is an exogenous variable independent of Pa(U,) and Zz,, j # The latent factors, as well as their possible regulatory relationships, are unknown.
In absence of interventions, the latent factors are sampled from the distribution PP and the mea- surements are derived through a decoder function g. i.e., u ~ Pu, + g(u). Interventions are assumed affect directly the latent variables U. rather than the observable x, and they can be either hard or soft (Pearl, 2009). In brief, hard intervention forcefully set the value of U; to a specific level, effectively severing any association between U, and its parents, while soft interventions solely modify the causal mechanism s,, altering the relationship between the variable and its regulators. Thus, under intervention 1, U is sampled from a new distribution w' Poz, while the decoder function remains unchanged, i.e., x² <- g(u').
In this context, the main goal of CRL is to identify a decoder function h and encoder function f such that h f(z) and fzz i. where ū "reconstructs" u as accurately as possible, while h approximates g. Additionally, one may be interested in identifying the regulatory (causal) mechanisms 8; among the U, components.
2

---

Published as a conference paper at ICLR 2025
Ahuja al. (2023) proved that be retrieved up to affine linear transformation, i..., it A u + с (Theorem 4.4 in Ahuja et al. (2023)). This requires two pivotal assumptions: (i) each intervention must targets single component of U, and (ii) the decoder function h is full rank polynomial. Notably, multiple interventions can target the same latent component U;, and, most importantly, the encoder function f is only required to be non-collapsing.
Zhang et al. (2024) further expand on this framework and reach the notable result that U can be retrieved up to permutation and scaling, i.e., 0, a, Uy + C (Theorem 2 in Zhang et al. (2024)). This result requires that the relationships U; + s((Pa(Ui) Zi) can be represented by Directed Acyclic Graph (DAG) with specific characteristics. and holds for both hard and soft interventions. Importantly. the authors proposed VAE-based architecture, the discrepancy- V.EE, that implements their theoretical results within a deep-learning framework. Describing the details of the discrepancy- VAE architecture is out of the scope of this work, however we note here few of its characteristics that are instrumental for the proposed SENA-discrepancy VAE:
The encoder f is implemented as a two-layer multilayer perceptron (MLP).
Once trained, the model provides two additional pieces of information: (i) a deep structural causal model (A, (0))²=) where the graph's 's adjacency encodes the parent of
each latent factor, while the matrix {s encodes the causal mechanisms (strength of interactions) (Pawlowski et al., 2020); and (ii) a map between each intervention and its target in the latent space, together with an estimate of the effect that the soft intervention has on
The variational nature of the discrepancy-Van VAE allows to predict the effect of unseen double perturbations, provided that each single perturbation is available during training.
The model trains both encoder f and decoder h only unperturbed cells, while the perturbed samples are solely used for deriving the effect of the perturbations on the deep structural causal model.
2.2 BIOLOGICAL. PROCESSES OR PATHWAYS.
biological process or pathway (BP) can be thought the of concerted biochemical reactions needed to perform a specific task within the cell (Kanehisa & Goto, 2000; Ashburner et al., 2000). In the context of this work, we loosely identify BP the genes contained within it, discarding information regarding other molecules or interactions. From this point of view, BPs can be simply thought as gene sets, where these gene sets can overlap or even contain one another.
2.3 CRL IN THE CONTEXT OF PERTURB-SEQ EXPERIMENTS.
In a Perturb-seg experiment, measurements x are single cell expression profiles, with each xy representing the expression of single gene i! Interventions are genetic perturbations in which one or multiple genes have their functionality inhibited (through, for example, a genetic knock-out, KO (Dixit et al., 2016)). In this sense. Perturb-seq perturbations represent hard interventions on genes: once knocked out, the level of functionality of the targeted genes does not depend anymore upon the other genes that usually regulate Two observations on genetic perturbations that are relevant for our main result:
Each perturbation likely affects several BPs at once. BPs are highly interconnected and genes are usually involved in several BPs at once.
Gene KOs (i.e., hard interventions) leads to soft interventions in BP activity. Biological systems are very resilient, partly due to high level of redundancy in their regulatory circuits (Reed et al., 2024). This means that following gene KO, other genes may partly assume the role of the suppressed gene, ensuring that the BP activity does not reach halt, even if it is somewhat impacted.
Herr we will assume that these values have been normalized and scaled to the point where they can be considered laying in R
3

---

Published as a conference paper at ICLR 2025
3 BIOLOGICALLY-DRIVER CAUSAL REPRESENTATION LEARNING
The CRL framework discussed in section 2.1 requires the decoder h be polynomial function. In contrast, much wider modeling flexibility is granted for the encoder f, which must simply be non-collapsing function. Thus, the encoder f can be built so it incorporates biological processes as prior knowledge.
To achieve this, we propose two-layer, masked multilayer perceptron (MLP) encoder, which we termed the SENA-ó (SparsE Network Activity) encoder (Figure 1). Let (BP1, BPx) be the gene sets corresponding to K BPs. Let org indicate the activity level of the k-th BP, summarizing to what extent genes within the corresponding BP are activated (i.e,, undergoing transcription). Then, the first layer of the proposed encoder connects the gene expression values x with BP activity levels o:
a = ( M)"
where W E R are the layer weights, a is the activation function, @ denotes element-wise multiplication, and M is mask matrix defined as:
{1 A
if gene i BPk, otherwise.
Mi.e
Each BP activity thus defined as a linear combination of the expression values of its respective genes. Unfortunately, it is known that the knowledge of the specific genes involved in BPs is seldom complete (Kunes et al., 2024). Thus, the tunable hyper-parameter À allows genes outside of the defined gene sets contribute to the BP activity if enough evidence of their involvement is present within the data. To this extent, À should be set to value small enough to discourage irrelevant contributions. Henceforth, we refer to this layer as the SENA layer.
The second layer follows VAE-type architecture where a fully connected linear layer with two heads (μ and o²) generates the exogenous variables Z, as:
2j ~ ~ where
B4 === ao
of =²=ad
Thus, the mean and standard deviation will be a linear combination of pathway activities a, weighted by the parameters d(m), 550), learned by the corresponding MLPs. Modeling each latent factor as a linear combination of BP activities, which we denote meta-pathway activities Zj. allows us to seamlessly combine the biological observation that each intervention may affects multiple BPs, with the CRL assumption (which provides identifiability guarantees) that each intervention must target only one latent factor (Figure 1). Modeling each latent factor as a single BP would set these two principles at odds with each other. We also note that each of the two layers of the SENA-ò encoder could be modeled as more generie, non-linear function, simply by adding intermediate layers. We opted for simpler architecture in order to prioritize interpretability over representational capabilities.
Most importantly, the SENA-ò encoder can be seamlessly plugged in the discrepancy-VA. VAE archi- tecture by substituting the original, fully connected MIP encoder. This modification guides the discrepancy- VAE VAAE architecture towards more interpretable subsets of the original solution space. We named this resulting model the SENA REeN--disrepancy -VAE (Figure 1).
Finally, we note that the SENA-6 encoder associates BPs activities to the Z,s, which in turn act as exogenous variables for the U.S, as noted in Section 2.1. The U,S are the actual causal latent factors involved in the causal graph, and the input for the polynomial decoder (see Fig. 1). To avoid confusion between these sets of latent variables, refer to the Z,s as meta-pathway activities (as each Z, incorporates the activity several pathways), and to the U,S as causal pathway archetypes. Appendix further elaborates the relationship between U,S and 2,8 through the causal graph A.
4

---

Published as a conference paper at ICLR 2025
SENNA-discrepancy--AE
discrepancy-VAE
SENA-6
Observational
Pathway activity
Pathway Archotype Activiyy
genes
Reconstructed Input
pathway Activity
Causal Pathway Archetype
Layer
2 puttee E-O" C-NO.1)
Polynomiall
Decoder
Pathway Mask
Causal Graph
MLP
Learnable
One Hot Encading
Meta-pathway
Intervened
Activity
archetype Activity
Pathway
Generated
Genes
Interventional
Interventional encoder
Decoder
Figure 1: Model overview. SENA disc-discrepancy- VAE modifies the encoder of discrepancy-VAn VAE to enforce biologically-drive training through a pathway-based mask M.
EXPERIMENTAL SETTINGS
To first assess the learning capabilities of our proposed architecture, we performed several ablation studies using the proposed encoder within simple autoencoder (AE) and variational-AF. (VAE) architectures. Our aim is to first assess whether the activity levels of the inferred (latent) BPs do indeed encode biological information (c.g., relevant activity changes are registered between perturbed and unperturbed data), and then to use these results to gauge the operational interval for X. We then compare the SENA disseepnccrepne---W VAA against the original discrepancy- VAE, over the task of predicting gene expression changes in novel combinations of perturbations, and finally analyze the interpretability of learnt latent causal factors.
4.1 DATA
We employ two large-scale Perturb-seq datasets, one collected on leukemia lymphoblast cells (K562 cell line) (Norman et al., 2019), termed the Norman2019 dataset, and a second one collected on acute myeloid leukemia cells (THP1 cell line) et al., 2022), termed the Wessels2023 dataset. The authors of Norman2019 targeted 12 genes known to affect the growth of the K562 cells, yielding total of 105 single-gene perturbation and 131 double-gene perturbations. The dataset underwent standard preprocessing steps for single cell data (filtering, normalization, and fog-transformation (Wolf et al., 2018)), leading total of 8,907 unperturbed cells (controls), 57.831 cells under the 105 single-gene perturbations, and 41.759 cells under the 131 double-geet gene perturbations. The Wessels2023 study underwent the same preprocessing steps, ending up include 424 unperturbed cells, 28 single-gen gene perturbations overall targeting 3036 cells, and 158 double- genee perturbations targeting 17592 cells. For all datasets we considered the 5,000 most variable genes in our analyses.
4.2 SELECTION OF BIOLOGICAL PROCESSES
Selecting an appropriate set of BPs is crucial, as the ideal selection should be sufficiently varied to include all active BPs in the system under study; and the same time, it is desirable to reduce the redundancy that usually characterize large sets of BPs. Following Ruiz-Arenas et al. (2024b), we considered the Gene Ontology (GO) BPs (Ashburner et al., 2000), and selected GO BPs with less than 30 genes. We then discarded those with more than half of their genes in common with other selected BPs, as well as those with low replicability (Ruiz-Arenal et al., 2024b). We further refine this selection by including only GO BPs that have least five genes represented in input dataset, and by removing those that are ancestors of other terms within our list. This multi-step selection ensures that the final BPs (mostly) non-overlapping and cover large variety of biological processes.
5

---

Published as a conference paper at ICLR 2025
4.3 IDENTIFYING ACTIVATED BIOLOGICAL PROCESSES
We exploit the architecture of the SENA-ò encoder to identify BPs that are activated under specific perturbations, In particular, we expect that the activated BP, should have its inferred activity level ag significantly altered with respect to unperturbed controls for those perturbations targeting genes within that BP. To measure this effect, we define the differential activation (DA) for BP; under intervention as DA = L o% where o and Of, are the values of activation function for BP; averaged over all perturbed and control cells, respectively. Because is a difference among mean values, t-test or similar inferential statistic can be used for assessing its statistical significance.
We then built two metrics for assessing to which extent the differentially activated BPs (i.e., statisti- cally significant DA values) are biologically meaningful. First, for each intervention p we define Wp as the set of BPs that contain the targeted gene W₂ = {BPx M, = 1}. The remaining (not affected) BPs are indicated as Wp. Intuitively, we would expect BPs containing the targeted gene to be the most affected by the intervention, while the other processes should only suffer indirect effects.
We then define the metric Hits@N, as the percentage of BPs in Wp that are ranked within the first N positions in terms of DAP- The parameter N is set 100 in analyses. Let Rig be the rank for BP; under perturbation p according to DA Then, Hits @N is defined as:
H
Σ kEW,
H S NJ.
(4)
W pl
Finally, we define the differential activation ratio (DAR) for a perturbation p as:
Wpl ΣkEW, DA DAR, W EREW, DA
(5)
This ratio contrasts the average activation for BPs directly affected by the perturbation against the average of the remaining BPs. Although computing this metric involves aggregating pathways with varying numbers of targeted genes, the imposed minimum of five genes per gene set and the definition of DA as an activation ratio make this metric potentially robust to intrinsic noise within the pathways. Note that both metrics require Wn and W, to contain at least one BP.
5 ABLATION STUDY
Due to the high sparsity infused in the SENA-5 encoder, one cannot assume that such encoder has good reconstruction capabilities while maintaining interpretability. Moreover, also seek to trade-off driven by the À parameter, Hence, in what
understand the wonnsstuccioonnnerrrtttbaee follows we assess the SENA-6 encoder by employing it in an AE and VAE architectures (with MLP as decoders in both cases), and compare it with a fully-connected encoder (denoted MLP) and two fi-regularized encoders with À as the regularization parameter. To perform fair evaluation, these architectures will each present two fully connected layers at the encoder, and the f encoders will only have the first layer regularized to imitate SENA's 's sparsity. To this end, we used the Norman2019 dataset. We evaluated the aforementioned architectures for several values of X: {0, 0.1, 0.01, 10-3). Overall, four different aspects were assessed: data reconstruction and generative capabilities (for VAEs), and interpretability and sparsity of latent dimensions, described next.
Data reconstruction and generative capabilities. We evaluated the reconstruction and generative capabilities of the proposed models by computing the test Mean Squared Error (MSE) and Kull- back Leibler divergence (DkL), respectively. The latter is only used the variational architectures.
Interpretability and sparsity of latent dimensions. We evaluated the interpretability by measuring how differentially activated (DA, see above) the affected neurons (i.e. BPs containing the knock-out gene) when compared the rest of the neurons after the SENA layer. Hence, compute the Hits@ 100 metric to measure the percentage of affected neurons the top 100 DA BPs. Moreover, we define the sparsity of a model the percentage of neurons with activation value (i.e., J Wi, smaller than 10 8, which measures the contribution of every input gene each factor after the SENA layer, and 7 refers to the mean expression across samples (cells). Reported metrics were computed on test samples.

---

Published as a conference paper at ICLR 2025
Results. Overall, enabling residual connections between genes and BPs in fully-connected fashion
(A = (100, 10 1103)) maintains paglooicallymmennnnge yielding reconstruction capabilities in par with the fully-connected MLP (Appendix IV Table 4 & Fig. 9-B). The reason could be that these models present an efficient use of the model weights (Appendix IV Fig. 8-B), underscoring the relevancy of gene-BP relationships in Perturb-seq data. Interestingly, higher values of lambda (A = 0.1) presented better reconstruction capabilities than the MLP (Appendix VII Table 4), at slightly sparser encoder (and hence, more interpretable) than the MLP (Appendix VII Fig. 8-B). On the other side, and as expected, À = 0 presented highly interpretable latent factors the cost of a significant drop in reconstruction capabilities. Additionally, our results show that the é--regularized MLPs does not perform well nor do they provide interpretable latent factors. Of note that when the analysis was perform only using the SENA layer, similar insights were obtained (Appendix VII Table 3 & Fig. 9-A). Finally, regarding the generative capabilities assessed on VAEs, the models based on SENA-ò encoder clearly outperforms other encoders on DKL. with lower values of X being the best performing ones (Appendix VII Tables 4 and 3, VAE-based column). These results were generated enforcing that every BP contains at least 5 genes. Different such thresholds modifies the number of considered BPs in the SENA layer (Appendix IV Fig. 10-B), slightly affecting the learning capabilities (Fig. 10-A) and performance time (Fig. 10-C).
latent factors (Appendix IV Fig. 8-A) while
6 LEARNING INTERPRETABLE LATENT CAUSAL FACTORS
We next benchmarked BENA-discrepaan VAE against its original counterpart, to assess the modeling and predictive capabilities of both models. This section focuses on the results obtained on the Norman2019 dataset, while Appendix V reports the results on the Wessel2023 data. For both datasets, we trained both models on the unperturbed and single-gene perturbations samples from Norman et al. (2019) (the latter are only used as a ground-truth -trth for the MMD loss). We also benchmarked GEARS (Roohani et al., 2024), a state-of-the-ait trtttrtertaaaaa approach for multigene perturbation prediction. Double-gene perturbations were set aside for evaluation purposes. We train both models across 3 different runs with default settings (Appendix F of Zhang et al. (2024)). Given the good results (in interpretability and reconstruction performance) obtained in the ablation study (Section 5), we varied the number of latent factors within (5, 10, 35, 70, 105), and the X for the SENA wiii-edepcereee---V varrepnnnyVVEE.e in {0, 0.13 (Appendix VII Fig. 12 shows gradients and mask (M) distribution across several À values).
6.1 PERFORMANCE BENCHMARKING
Table reports the results of the comparison, where MMD (Max Mean Discrepancy (Gretton et al., 2012)) measures the difference between the generated and true double-perturbation distributions. We report the average MMD over all 131 double-gent geuble-gene perturbations. Additionally, MSE indicates the reconstruction error for control samples during training, DKL is the variational loss (Kingma & Welling, 2014), and LI == ||A|] represents the sparsity of the deep structural causal model.
Table I: Benchmarking SENA-discrepancy -VAE and discrepancy- VAE on double perturbations prediction. Values reported as mean = variance computed on 5 runs with different initializations.
Latent Dimension
Encoder
Metric
105
70
35
10
MMD.
1.59811
0012110
1.73486
0.022115
1.98993
00030033
2.43440
0.000570
2.53237 + 0.077662
Original MLP
MSEL
0.02152
0.00155
0,02298
0.000554
0.02499-
0.0000s
0.02699
:::
0.02792 1 0000015
KLDJ
0.00022
+ 0.00005
0.00021
+ 0.000055
0.00021
+ 0,00000
0.00024
+ 0.002033
0.00030 + 0.000027
L14
0.06097
0003718
0.06934
0000055
0.06714
0,0066s
0.07201
0.009514
0.08319 + W000090
MMD
1.58489
0.006100
1.75332
0,004884
1.94984
0.027583
2.49881
0.066160
2.60439 1: 0.100705
SENA-6X-0.1
MSEL
0,02134
+ 0000047
0,02300
00000120
0.02462
0.000041
0.02688
+ 000000066
0.02805 + MAM00135
KLDE
0.00019
0000001
0.00019
0.00006
0.00019
0.0000n
0.00020
00000007
0.00020 + 00000522
LL
0,05243
00015441
0.05323
0.07308
0,001523
0.08589
0,06856 +
MMD 1.74588
00038799
1.90425
2.22812
0.154388
2.62393
2.90044 + 174605
SENA-82-0
MSEL
0.02312
+ 0.000488
0.02460
+ 0.000088
0.02634
1: 0,000066
0.02800
1 0.000236
0.02879 1 0.001447
KLD
0.00018
0000002
0.00018
00000000
0.00019
6.0000m
0.000221
0.002017
0.00022 + 0000021
L11
0.05415 18 0.001903
0,05749
0.000457
0,06730
0,000 195
0.08651
+ 0.025800
0.09486
1: 0.08236
GEARS
MMD
14.9420 + 233857
12.6036
0.745302
13.2590
0.559124
12.9774
1: 0.751373
15.3099 1: 0.191833
7

---

Published as a conference paper at ICLR 2025
Interestingly, and despite the restrictions imposed by the SENA-ò encoder that could potentially decrease the SENA -iscrrpaacc--.... VAE representational capabilities, the proposed model outperformed the MLP encoder for some latent dimensions in terms of MSE and MMD computed on unseen double perturbations for small values of À (0.1). Moreover, setting À = 0 allowed the JENA-diicrepancy- VAE to surpass the original MI.P encoder on the DKL. metric, while the optimal model for causal graph sparsity (LI) varied with latent dimensions. These results, which align with those of the ablation studies, highlight the potential of SENA-discrepancy-VV failed to properly model the evaluated double perturbations. Note however that GEARS does not provide a causal graph nor a generative model (details in Appendix VI).
On the other hand, GEARS
6.2 VISUALIZING SENA DISCREAACYE VAE LATENT FACTORS
We first investigated the association between perturbations and latent factor activation (Fig. 14). Both models tend to activate few latent factors. Specifically, the discrepancy- VAAE model activate 8 to 9 factors across all perturbations when 35 or more latent factors are included in the model. These numbers decrease to 6 and 4 when 10 and 5 latent factors are available, respectively. At the same time, more than half of the perturbations are assigned to only 1 or 2 latent factors, creating a quite unbalanced mapping. The SENA-discrepaan that relatively few latent factors are needed for capturing the changes induced by perturbations, while the remaining latent factors assist representing the overall distribution of gene expression data.
VAE follows a similar pattern. This seems to indicate
Interpretation of the SENA-discrepancy, VAE latent factors. The proposed model offers the possibility of inspecting its en- coder for deriving the BPs composing the latent factors. By construction, each per- turbation will target a single latent factor U,. which enable us to associate each BP to the intervention with the largest differential activation value. Only significant differen- tial activation values are taken into account (ranked within the top 1 1 in absolute value, and false discovery rate (FDR) S 0.05 via two-tailed t-test with BH correction). Figure 2 represents the causal graph for the latent factors associated to at least one BP for the SENA -discrppaccc-aaii VAAE model with 105 la- tent dimensions and À = 0. Ten edges with the highest coefficients in absolute value are reported for readability. Latent factors are represented
65
69
15
53
Figure 2: SENA SENA-discrepancy- VAE causal graph on the Norman2019 data. Latent factors represented by word clouds of associated BPs. Arrows indicate causal influences.
word cloud of BPs (i.e., a graphical representation of the terms fre- quencies within the BPs names). Table 7 (Appendix VII) reports the number of perturbations and BPs assigned to each factor and Table 9-10 (Appendix VII) shows the mapping between BPs and selected latent factors. It worth noting that the inferred causal graph is robust across latent dimensions and As, where most of the inferred connection are maintained (Appendix II).
Latent factor 15 is targeted by perturbations on the JUN gene, and is associated with the activity level of GO:0050665, "hydrogen peroxide biosynthetic process". While JUNE is not included in GO:0050665, this gene is known react with over-expression to oxidative stress (Vandenbroucke et al., 2008). For latent factor 69, the targeted PTPN13 gene is known to be involved in several tumors (Mcheik et al., 2020), thus it is not surprising to found it associated with BP related to blood vessel formation. Interpretation of other factors requires careful inspection of their associated BPs (see Appendix VII Tables 9 & 10). For example, most of the BPs in latent factor 65 are associated to tissue development, while latent factor 12 contains several BPs related to protein activity.
Upon inspecting the connections on the causal graph, a first important connection is the one between factor 15, "hydrogen peroxide biosynthetic process" process""Appeenn "AAppendix VII. Table 9, third latent factor) which causes factor 69, "endothelial cell morphogenesis" (Appendix VII. Table 10, last row). It is well known that hydrogen peroxide stimulates endothelial cell proliferation (Stone & Collins,
8

---

Published as a conference paper at ICLR 2025
2002; Anasooya Shaji et al., 2019). Thus, causal graph captured this regulatory relationship in a fully unsupervised, data-driven way. In turn, factor 53 causally influences factor 15, and factor 53 contains the biological process "catechol-containing compound biosynthetic process" (Appendix VII, Table 10, second element of factor 53). It is well known that H202 can be produced by the metabolism of catecholamines (Noble et al., 1994; Seregi et al., 1982). An connection exists between latent factor 69 and latent factor 2, with the latter including "negative regulation of endothelial cell apoptotic process" (Appendix VII, Table 9, second row) among its biological processes. Taken together, these findings provide evidences for the correctness of our approach and its capability of recapitulating known biological causal relationships.
more direct
6.3 SENN-DISCREPANCY- VAE CAPTURES BIOLOGICALLY MEANINGFUL PATTERNS
We next evaluate the ability of the proposed encoder to maintain hologically-drivve factors (see Appendix I). For this evaluation, we focus on the Norman2019 dataset, and we set the latent space dimension to 105, one for each single-gene perturbation in the dataset. We then evaluated the 37 knocked out genes that were present at the input. For each of these, we computed the DA score across all BPs (after the SENA layer) and found that those including the targeted gene reported higher DA on average (Appendix VII Fig. 13). Since each perturbation presented a different number of affected BPs, next focused the perturbations with the largest amount of targeted BPs (i.e., 7), and evaluated the significance (statsannotation package (Charlier et al., 2022)) of the DA among affected and not U test with BH p-value correction). Fig. 3-A shows the aforementioned analysis for the knock out genes LHXI, SPII and TBX3. This analysis highlighted the perturbation samples from LHXI and TXB3 KOs those presenting highly differentially activated BPs, validating the capacity of the proposed encoder to identify biologically menninggf factors.
affected BPs
B
DAR
13
10"
10"
17"
10
1900)
Knockout
Lateni Dimansion
Layer
Figure 3: SENA, l and o² layers interpretability analysis. A. DA score for the three perturbations presenting the highest number of affected BPs among W. and W.- B. DAR of the 37 analyzed perturbations when varying the latent space dimensions, For every dimension, genes with the highest DAR are highlighted. C. DA score across the evaluated perturbations at the output of the and o² layers among affected and not affected BPs. ns: non-significant, ****; p-value 10-4
To analyze the robustness of the SENA- dire--ddsseeepcy--yye VAE's encoder, we repeated the above analysis across several latent space dimensions. computing the DAR (Eq. 5) between those BPs containing the targeted gene and the rest. Once again, we found that almost every evaluated knock out gene reported DAR 1 along evaluated latent dimensions (Fig. 3-B). It is worth highlighting TMSB4X and LYL.I, 10 consistently, indicating that the mean difference in activation of the SENA -layrr neurons among perturbation and control samples for the BPs containing targeted genes was > 10x times greater than for the rest. This underscores the capacity of SENA SENNA-disrrepane--Aal factors.
which reported a DAR >>
drive the training process while maintaining
Interpretable μ and x² layers. We next assessed how the above shown interpretability is propagated through the SENA -ENA--iicrepaacc--- VAE encoder. To this end, we performed the DA analysis on the output of SENA-ò - encoder (Figure 1), where we measured the contribution of affected and not affected BPs from the previous layer every neuron j the fu and o² layers. To this end, we define the DA
score for perturbation and BP at the j-th neuron of and o² layers as (DA?), ) = |8kz baall E 0% where superscripts p and e denote perturbed and controlled activities, respectively.

---

Published as a conference paper at ICLR 2025
Fig. 3-C shows the DA score on μ and o² layers for affected and not affected BPs, following the same significance tests performed above. Again,
A
Interventions
GO:0038065
this highlights that the MeeA--disreppncy-.. VAF the
CEBPA
encoder maintains biologically-mannine
UMMP2
COL1A1
GLB1L2
SENA-8 encoder (Fig. 1), yielding interpretable exogenous variables Zi. which we denoted the meta-pane pett-ppttwway activities and also contributes to this differential activation. We next analyzed the biological significance of the meta-pathway activ-
2
TEXS
-2-
-15
TP73
-2
-2
UMAP1
O
2
-8
-2
C
D
ities, since these connections are learned in a data- driven manner. For this, we measured, through a permutation test, the contribution of the level 2 GO pathways (i.c., the parental GO terms of the used BPs in the SENA layer) to each of meta- pathway activity. Interestingly, multiple meta pathways were significantly associated to few level 2 pathways, underscoring our model capabil- ities to learn hiologically-mennnngfff meaningful patterns at both high (BPs) and broad (meta-pa aneta-pathway) gran- ularities (Appendix III.)
GO-0006633
G0:0070141
2-
0-
0-
-05
42
-12
-10
-30
-1
-25
2
E
F
60:0010714
GO:0008595
4-
80
0-
0-
Finally, we performed differential activation score analysis on the Norman2019 dataset after training SENA distr-pdcceepacy---e VAE. We selected the top 6 largest DA scores, which belonged to unique KO genes and gene sets, respectively (Table 8). Fig. 4-A shows the UMAP components of all intervened cells (in the input gene space) across the aforementioned genes, while Fig. 4 B-F depicts those cells colored by the DA score that each cell has the evaluated gene set (GO term). Surprisingly, from the top 6 DA scores, we found that only GO:0038065 is initially targeted by COLIA1, while the remaining gene sets are reporting specific-highlee
-10
-20
-25
-6
Figure 4: DA score of 6 most significant (KO, BP) pairs. A. UMAP of all intervened cells across the KO genes presenting the 6 largest DA scores. B-F. UMAP from A colored by the differ- ential activation score on each cell across genes within each evaluated BP.
on the cells belonging to the intervened KO without being directly targeted. For instance, the gene set GO:0006833 (Fig. 4-C) is activated by the TBX3 gene (same with GO:00107 14 and CEBPA) without being explicitly encoded in the SENA-6 encoder. These underscore the potential of SENA-discrepancy- -VAE VAE to naturally learning biologically-drive patterns without specifically enforcing them.
DISCUSSION AND CONCLUSIONS
In this work we have demonstrated how biological processes can be used as prior knowledge in the context of causal representation learning. The resulting model, SENA-din discrepancy- VAE", is on par, or even outperforming in specific scenarios, in terms of predictive capabilities with the original discrepancy- VAF, while at the same time producing embeddings that can be easily inspected for assessing their biological meaning.
Among the the several findings reported in this study, it is striking that both models tend to assign most interventions to a small number of latent factors (see Fig. 14 for perturbation-to- associations in the Norman2019 data). Reasoning in terms of biological processes helps understanding why. The theory behind discrepancy- VAE requires that each intervention must be assigned to single factor. Thus, this factor must represent all BPs affected by that intervention. If two (or more) interventions affect overlapping sets of BPs, then by necessity all these BPs must be mapped to the same factor. Overall, it may be argued that assuming that each intervention targets single latent factor does not allow CRL methods to thoroughly disentangle the interplay between BPs, perturbations and latent factors. Thus, future CRL works should attempt overcome this assumption.
factor
Python package, including data and code for reproducibility: github.com/Mi LBM--LLLa//EEAAA
10

---

Published as a conference paper at ICLR 2025
REFERENCES
Kartik Abuja, Divyat Mahajan, Yixin Wang. and Yoshua Bengio. Interventional causal representation learning. International conference on machine learning. pp. 372 407. PMLR, 2023.
Chinchusha Anasooya Shaji, Bobby D) Robinson, Antonia Yeager, Madhava R Beeram, Matthew Davis, Claire I. Ishell, Jason H Huang, and Binu Tharakan. The tri-phasic role of hydrogen peroxide in blood-bre brain barrier endothelial cells. Scientific reports, 9(1):133, 2019.
Michael Ashburner, Catherine A Ball, Judith A Blake, David Botstein, Heather Butler, J Michael Cherry, Allan Davis, Kara Dolinski, Selina Dwight, Janan Eppig, et al. Gene ontology: tool for the unification of biology. Nature genetics, 25(1):25-29) 2000.
Florian Charlier, Mare Weber, Dariusz Izak, Emerson Harkin, Marcin Magnus, Joseph Lalli, Louison Fresmais, Matt Chan, Nikolay Markov, Oren Amsalem, Sehastian Proost, Agamemnon Krasoulis, getzze, and Stefan Repplinger. Statannotations, October 2022. URL https: //doi http//iii.o////////////////// org/10.
5281/zenodo. 7213391.
Yongin Choi, Ruoxin Li, and Gerald Quon. sivae: interpretable deep generative models for single-cell transcriptomes. Genome biology, 24(1):29, 2023.
Haotian Cui, Chloe Wang, Hassaan Maan, Kuan Pang, Fengning Luo, Nan Duan, and Bo Wang. scgpt: toward building a foundation model for single-cell multi-omics using generative ai. Nature Methods, pp. 1-11, 2024.
Atray Dixit, Oren Parnas, Biyu Li, Jenny Chen, Charles Fulco, Livnat Jerby-Arnon, Nemanja D Marjanovic, Danielle Dionne, Tyler Burks, Raktima Raychowdhury, et al. Perturb-seq: dissecting molecular circuits with scalable single-cell raa profiling of pooled genetic screens. cell, 167(7): 1853-1866, 2016.
Thomas Gaudelet, Alice Del Vecchio, Eli M Carrami, Juliana Cudini, Ghantriolnt-Andren Kapourani, Caroline Uhler, and Lindsay Edwards. Season combinatorial intervention predictions with salt & peper. arXiv preprint arXiv:2404. 16907, 2024.
Luke Gilbert, Max A Horlbeck, Britt Adamson, Jacqueline Villalta, Yuwen Chen, Evan H Whitehead, Carla Guimaraes, Barbara Panning, Hidde I. Ploegh, Michael C Bassiik et al. Genome- scale crispr-mediated control of gene repression and activation. Cell, 159(3))647-661, 2014.
Arthur Gretton, Karsten M. Borgwardt, Malte J. Rasch, Bernhard Schölkopf, and Alexander Smola. A kernel sao-la test. Journal of Machine Learning Research, 13(25))723-73,, 2012. URL http:// jmlr. org/ papers /v1 /tt3////////////eee ggwgttttttnaaaa 2a. html.
Gilles Gut, Stefan Stark, Gunnar Rätsch, and Natalie Davidson. Pmvae: Learning interpretable single-cell representations with pathway modules. bioRxiv, pp. 2021-01, 2021.
Minoru Kanehisa and Susumu Goto. Kegg: kyoto encyclopedia of genes and genomes. Nucleic acids research, 28-1): 2000.
Diederik P. Kingma and Max Welling. Auto-Encoding Variational Bayes. In 2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings, 2014.
Russell Z. Kunes, Thomas Walle, Max Land, Tal Nawy, and Dana Pe'er. Supervised discovery of interpretable gene programs from single-cell data. Nature Biotechnology, 1089) Jul 2024. ISSN 1546 1696. doi:
URL httpst /doi http://ooooo...../////o/// org/10.
023 - 01940-3.
Romain Lopez, Jeffrey Regier, Michael B Cole, Michael Jordan, and Nir Yosef. Deep generative modeling for single-cell transcriptomics. Nature Methods, 15(12):1053-1058 2018.
Romain Lopez, Natasa Tagasovska, Stephen Ra, Kyunghyun Cho, Jonathan Pritchard, and Aviv Regev. Learning causal representations of single cells via sparse mechanism shift modeling. In Conference Causal Learning and Reasoning, pp. 662-691. PMLR, 2023.

---

Published as a conference paper at ICLR 2025
Mohammad Lotfollahi, Felix A Wolf, and Fabian J Theis. scGen predicts single-cell perturbation responses. Nature Methods, 16(8):715-721, 2019.
Mohammad Lotfollahi, Sergei Rybakov, Karin Hrovatin, Soroor Hediych-zadeh. Carlos Talavera- López, Alexander V. Misharin, and Fabian J. Theis. Biologically informed deep learning to query gene programs in single-cell atlases. Nature Cell Biology, 25(2):337-350, Feb 2023. ISSN 1476-4679, doi: 10. 15538/44155--02-000772. URL https: /tw./iiooooooooooooo wtt.:.........oooo 10388/ s 2-1156-0222-00722
Soha Mcheik, Leticia Aptecar, Peter Coopman, Véronique D' D'Hondt and Gilles Freiss. Dual role of the ptpn tyrosine phosphatase in cancer. Biomolecules, 10(12), 2020. ISSN 2218-273X. doi: 10.3390/bim101126659) URL https: /twp//dwpppp.....m////// com/ 2218 273X/10/12/1653
Seyednami Niyakan, Xihaier Luo, Byung-Jun Yoon, and Xiaoning Qian. Biologically interpretable vae with supervision for transcriptomics data under ordinal perturbations. bioRriv, pp. 2024-03, 2024.
Paul G Noble, Jack P Antel, and Voon Wee Yong. Astrocytes and catalase prevent the toxicity of catecholamines oligodendrocytes. Brain research, 633(1-2):83-90) 1994.
Thomas M Norman, Max A Horlbeck, Joseph M Replogle, Alex Y Ge, Albert Xu, Marco Jost, Luke A Gilbert, and Jonathan S Weissman. Exploring genetic interaction manifolds constructed from rich single-cell phenotypes. Science, 365(6455);786-799) 2019.
Nick Pawlowski. Daniel Castro, and Ben Glocker. Deep structural causal models for tractable counterfactual inference. In Advances in Neural Information Processing Systems 33 (NeurIPS 2020), 2020. doi: (0.555/3495724.449559999 URL https: /dl. acm orq doi//0..11000/ 5555/ 3495724. 3495797.
Judea Pearl. Causal inference in statistics: An overview. Statistics Surveys, 2009.
J. Michael Reed, Benjamin E. Wolfe, and L. Michael Romero. Is resilience a unifying concept for the biological sciences? iScience, 2764)80 2024. ISSN 2589-0042. doi: https://doi. org/10. 101..0000/i...ii.200....... 2024. 109478. URL https: www..wwwwoooooom sci encedi rect com/ science/ (TT.cc/e//////233000240099 0002222889904220099990
Yusuf Roohani, Kexin Huang, and Jure Leskovec. Predicting transcriptional outcomes of novel multigene perturbations with gears. Nature Biotechnology, 42(6):927-935, 2024.
Carlos Ruiz-Arenas, Irene Marín-Goñi, Liewei Wang, Idoia Ochoa, Luis A Pérez-Jurado, and Mikel Hernaez. NetActivity enhances transcriptional signals by combining gene expression into robust gene set activity scores through interpretable autoencoders. Nucleic Acids Res., 52(9):e44, May 2024a.
Carlos Ruiz-Arenas, Irene Marín- Monin Liewei Wang, Idoia Ochoa, Luis A Pérez-Jurado, and Mikel Hernaez. Netactivity enhances transcriptional signals by combining gene expression into robust gene set activity scores through interpretable autoencoders. Nucleic Acids Research, 52(9):e44 e44, 2024b.
Bernhard Schölkopf, Francesco Locatello, Stefan Bauer, Nan Rosemary Ke, Nal Kalchbrenner, Anirudh Goyal, and Yoshua Bengio. Toward causal representation learning. Proceedings of the IEEE, 109(5):612) (1-5)6612-634 2021.
Lucas Seninge, Ioannis Anastopoulos, Hongxu Ding, and Joshua Stuart. Vega is interpretable generative model for inferring biological network activity in single-cell transcriptomics. Nature communications, 12(1):5684, 2021.
András Seregi, Péter Serfözö, Zsuzsanna Mergl, and András Schaefer. On the mechanism of the involvement of monoamine oxidase in catecholamine- -tamiieettimultted prostaglandin biosynthesis in particulate fraction of rat brain homogenates: role of hydrogen peroxide. Journal of Neurochemistry, 38(1):20-27, 1982.
James R Stone and Tucker Collins. The role of hydrogen peroxide in endothelial proliferative responses. Endothelium, 94)8231- 2002.
12

---

Published as a conference paper at ICLR 2025
Alejandro Tejada-Lapuerta, Paul Bertin, Stefan Bauer, Hananeh Aliee. Yoshua Bengio, and Fabian J. Theis. Causal machine learning for single-cell genomics, 2023. URL httpttt/aaxx.i..oom
0ts/2310..4435,
Korneel Vandenbroucke, Steven Robbens, Klaas Vandepoele, Dirk Inzé, Yves Van de Peer, and Frank Van Breusegem. Hydrogen peroxide-induced gene expression across kingdoms: a comparative analysis. Mol. Biol. Evol., 25(3):007-516, March 2008.
Hans Hermann Wessels, Alejandro Méndez-Mancilla, Yuhan Hao, Efthymia Papalexi, William M. Mauck, Lu Lu, John A. Morris, Eleni P. Mimitou, Peter Smibert, Neville E. Sanjana, and Rahul Satija. Efficient combinatorial targeting of ma transcripts in single cells with cas ma perturb-seq. Nature methods, 20:86. 2022. ISSN 15487105. attps://pmet ncbi nlm. nih. httphh/wwwarriceeeee//0000330000
doi: 10.1038//41192-022-00777
URL
F. Alexander Wolf, Philipp Angerer, and Fabian J. Theis. Scanpy: large-scale single-cell gene expression data analysis. Genome Biology, 19(1), February 2018. ISSN 1474-760X. doi: 10. 1186/
$13059-017-1382-0 URL http::/tpeee /xx doi.
Jiaqi Zhang, Kristjan Greenewald, Chandler Squires, Akash Srivastava, Karthikeyan Shanmugam, and Caroline Uhler. Identifiability guarantees for causal disentanglement from soft interventions. Advances in Neural Information Processing Systems, 36, 2024.
Yifan Zhao, Huiyu Cai, Zuobai Zhang, Jian Tang, and Yue Li. Learning interpretable cellular and gene signature embeddings from single-cell transcriptomic data. Nature communications, 12(1): 5261, 2021.
13

---

Published as a conference paper at ICLR 2025
APPENDIX
I INTERPRETABILITY OF LATENT FACTORS AND CAUSAL GRAPH THROUGH OUR PROPOSED SPARSE LAYER
In the variational autoencoder proposed at Zhang et al. (2024), the exogenous variable Zj is sampled from normal distribution, where the mean and standard deviation of this distribution is defined by the fully connected layers in their encoder. In the proposed SENA -iin---aerppancy--. Serepancy-VVEEE (Fig. 1), the mean and standard deviation will be a linear combination of pathway activities a. weighted by the parameters ) 0 learned by the corresponding MLPs. Thus,
My = a a of a² o
(6)
which then define the meta-pathway activities 2j as 2j Ni1;; of). This allows the expectation of Zj to be interpretable, as
B(Z,) = 1 = 8'
However, we would also like this interpretablity to hold when going from the exogenous variables (the meta-pathway activities) to the causal factors (causal pathway archetypes). The latter are
defined as U z (I A)- where
L a (1 A)
Σ 1=0
AK),
(8)
according to the Neumann series, and given that represents the adjacency matrix of Direct Acyelic Graph, with L being the largest path (hence, AA = 0, k = K + .). Here A defines the causal relationships in the latent space. Therefore, the j-th causal factor can be expressed
U==
where L5 is the jth column of L, and encodes the number of path of at most length K that ends at node j in the causal graph. Hence, the expectation the causal factor U, is given by
RR(U,) =R(Z* L2)
=F(Z") Lj
.e. Lj
(9)
=a". A Lj
(10)
= =a= -5/40
where Eq. (9) and Eq. (10) from Eq. (7), and A (H) is the (learned) linear mapping between the pathway activities a and the meta-pathway activities Z. Thus, 5(th) (linearly) maps the causal latent factors with the pathway activity scores through the learnt causal structure L.j. providing the mechanism for the interpretation of the latent factors, termed in our work as the pathway archetype activities.
Finally, we experimentally validated, using the Norman2019 dataset, that Eq. (9) holds for both the original discrepancy- (MLP), and the proposed SENA-ò model for both x= (0, 0.1). We used all available unperturbed cells (ctrl) and 9 randomly- chosenn perturbed cells, setting the latent dimension to the number of available perturbations in the dataset, i.c. 105. To this end, first computed B(U;) for every latent dimension j E {I...., 1053 by forwarding the cells and averaging (over 10,000 realization of Z NN o²)) the obtained pathway archetype scores (i.e., the causal latent factors Us). On the other hand, we multiplied the BP activity scores (a) with the BP-to- meta- neta-pathway mapper () and the causal mechanism (L. Eq. 8). Both computations should be equal according Eq. (10). Fig. 5 depicts for every type of perturbation and model (MLP, SENA-ôX20 and SENA-6x both terms of the equality, as well as the Pearson's 's correlation. There a perfect correlation among these terms, and this patterns is maintained across models and perturbations.
14

---

Published as a conference paper at ICLR 2025
2- 100
Tree TAN
TAM
TAM
1/0
TAM
LAM
TUM
1,00
-3
-25
F(U)
11100
11/0
10
14100
1100
25
1.000
10
1100
14100
14100
11100
+100
15
-20
10
-25
20
-05
12
-10
-5
+1.
-10
-35
-15
-15
1.0
1.00
11.00
2.
35-
-12
50
50
50
-0.5
-35
-88
-15
-25
-86
Figure 5: Experimental validation of derived Eq. 9. A-C. Experimental analysis of derived Eqs. and 10 for the original discrepancy-VAR. (MLP) architecture (A), SENA-ôx=0.1 (B) and SENA-ôx=0 (C), respectively. Here x-axis represent the expected latent factors extracted by forwarding cells through the trained model for each selected cell type and y-axis represent the (linear) mapping between BP activation scores and pathway archetypes (i..e., latent causal factors).
15

---

Published as a conference paper at ICLR 2025
II STUDY ON INFERRED CAUSAL GRAPH ROBUSTNESS
We evaluated how stable are the (directed) edges from the inferred SENA SENNA-iisrepaancy- Deaiiiscrrpaa--VVVEEE. VAE causal graph across various A's and latent dimensions(Fig. 2). To this end. inferred the causal graph for À = {0, 0.1, 10-2, 10-3} and latent dimensions {5, 10, 35, 70, 105} using the Norman2019 9 dataset. We then computed the edge consistency for each graph as the ratio of the frequencies of the most and least frequent sign across the evaluated A's. For instance, if an edge has been consistently positive across s, it would present an edge consistency of (e.g., Fig. 6-B, U12).
When analyzing the edge consistency of the inferred causal graphs (e.g.. Fig. 6 A-B, depicted for latent dimension 5), most edges had a consistency above or equal to 75% while edge weights close to presented low consistency across X's. Interestingly, across the different tested hyperparameters, most edges had a coefficient of variation (CV) 2 (Fig. 6 C-F), indicating perfect edge consistency, Importantly, the large majority of edges exhibited high confidence (CV 2 D)(Fi, 6-H).
A
E
G
20-
6,22
Wioe Volue
11
-0.30
B
H
Ledee
163
80
0.
302.5%
Caal (10000100)
Figure 6: Edge robustness analysis for the inferred causal graph. A. Edge values (mean of the inferred graph upper triangular. Results averaged across {0, 0.1, 10- 10-3). B. Edge consistency across the aforementioned X's. C-G. Coefficient of Variation (mean over standard deviation) histogram for every edge of the upper triangular matrix of the inferred causal graph across several runs using different latent dimensions ((5, 10, 35, 70, 105 H. Cumulative distribution of the edge consistency.
16

---

Published as a conference paper at ICLR 2025
III STUDY ON HIGH-LEVEL-ORDER AGGREGATION OF BIOLOGICAL PATHWAYS
This section analyzes how the connections between the BPs activity after the SENA layer and the exogenous variables (meta-pathways) that then define the causal latent factors (pathway archetypes) are potentially biologically-drin end, we relied on the level 2 Gene Ontology pathways, that encompass most of the used lower- level BPs in our model to evaluate whether the meta- secta-patwway activity encodes high-level biological processes. For the sake clarity, for this section we would refer to the aforementioned level pathways as I.2BPs and keep using BPs to refer to the ones we model at the output of the first layer SENA -iseeeeierypaacc-- VAE.
despite the weights being learned in a data-driven way. To this
We selected all L2BPs containing at least 10 BPs (from our set of filtered BPs, i.e. 454 if not stated otherwise), yielding a total L2BPs with different sizes (Table 2). In order to measure the association between the inferred meta-pathways and the L2BPs, we computed the contribution that each L2BP has on each meta- peeta- score and perform a permutation test to evaluate whether the association is significant. We now describe the process for computing these contributions. Be 2j the activation score of the meta-pathway j. this activation score is given by the linear combination of previous layer's activation (BPs' activation scores) and the learned weight matrix A of that layer.
Lets define L2BP, l as the set of BPs within the k-th L2BP. Since the activation score for every meta-pathway j can be expressed
2j
Σ REL2BP&
BP; 8ij
Σ
4k E {1, (1,..., [L2BP]] BP. * 0,y
I' #L22BP
we can compute the contribution of the k-th HIP to the j-th latent factor as:
EiELAAP BP:
Ckj
2j
To measure if this contribution is significant, for a given I.2BP and a given meta-pathway, we permuted 1000 times the BPs (maintaining the size) and computed the Mann-Whitne one-sided test to obtain the p-value between permuted and true BPs associated to the L2BP, correcting by the number of tests performed (Bonferroni correction). These allowed statistically significant contribution of the specific L2BP (and associated BPs) to the activation score of the given meta- pateaay. We performed this analysis on SENA-6 for {0, 0.1} and, for the sake of simplicity, we set the number of latent dimensions to 35.
measure there is
Fig. 7 A-B shows the histogram of permuted vs true contributions for every L,BBP on the first meta-patinn meta-pathway factor and Fig. 7-C depict the distribution of corrected p-values for ev- cry meta-pathway and L2BP, where blanks represent non-significant contributions (corrected p-value < 0.05). Interestingly, there is significant contribution for every L2BP across sev- eral meta-pathway factors, which may indicate that there are potential clusters of them encod- ing true high-level biological processes. Also, these results highly similar across the evalu- ated À values. Note that this analysis did not discriminate across perturbation types, hence all cells were forwarded (and averaged) through SENA's model compute the activation scores.
Table 2: Aggregation according to the Level-2 Biological processes of the Gene Ontology structure.
Level GO Term
#Genesets within
Description
GO:0022414
13
Reproductive process
GO:0002376
21
Immune system process
GO:0051179
Localization
G0:0032501
39
Multicellular organismal process
G0:0050896
41
Response to stimulus
GO:0008
Metabolic process
GO:0032502
64
Developmental process
60:001987
141
Cellular process
GO:0065007
Biological regulation
17

---

Published as a conference paper at ICLR 2025
CO: SCOSER H4)
USSC087
CO 02011790
CO 2522907
0085627
GC 0332500 (84)
DD: 2222414 (13)
CC 0002379 1211
CD 0006152 (83)
1.00 (0/40)
1:0000
1000
Score
2omm
Sooe
(141)
200
1000--00
30
10030-0000000 (4)
X.03
Score
Some
Don
D- 2:000000
2*2750025 (!T))
pa 1030x 00 1N
pa 20046.7 (TT)
50m
500
220
A 0.0
0.1
GO NOORAOR (P1)
GO 000000110 (93)
(141)
50
Goma 00 dhg
GOOKERSDI
32
GOUNERSOO (04)
DOUNSSERE 211)
Figure 7: Study on genesets aggregation at the latent factor level. A-B. Permutation test on the latent factor contribution for level 2 genesets and À = 0.1 (B). Results are shown for the first latent factor C. Heatmap depicting the corrected p-value for every level geneset and latent factor SENA-6 when
random aggregations of genesets, for À (A)
(left) and 0.1 (right). Number of BPs inside every level 2 GO term are shown in brackets next to the terms' name.
18

---

Published as a conference paper at ICLR 2025
IV
ABLATION STUDIES ON NORMAN2019'S DATASET
A
B
0.7
SENA.AE
SENA-G-AE
0.6
SENA-6-VAE
0.7
1 Spprsiiyy
0.5
H
0.4
0.5
0.4
0.3
0.3
0.2
0.2
0.1
0.1
10
10
10
0
SENA (w-0)
L1 SENA SENA SENA (0-1000-11000000000001
À
Figure 8: Ablation studies interpretability and sparsity. A. Percentage of affected gene sets in the top 100 DA BPs for several SENA-based architectures and values. B. Sparsity evaluation according to Eq.??
-B-M4P (1
Test (x100)
Test MSE (x10")
50
160
150
200
250
50
100
150
200
250
Epoch
Figure 9: Ablation studies for AE-type architecture. Test MSE evaluation for AE-based architec- tures for SENA (A) and SENA-ò (B) encoders.
19

---

Published as a conference paper at ICLR 2025
A
B
169
Mirku/een
réase
Acrounn
10
Tes1
10
Mixirum number of gunis witton each garseett
10"
50 75
100
125
150
175
200
Epoch
Minimum rumber genes with ench gennent
Figure 10: Ablation studies of the number of BPs used in SENA-ôX=o- On each experiment, minimum number of genes per BP is imposed, which reduces the total number of BPs used after the SENA layer. A. Test KLD as function of the minimum number of genes within eachBP. Summary of the number of BPs and the gene-to-BP connections. C. Training time-per-epoch results. Results are averaged across different seeds.
Table 3: AE and VAE-based evaluation for SENA across seeds. Methods are sorted by sparsity.
AF-based
VAF-based
Method
Test MSE (x 10
Test MSE (x10 3)
Test DgL. (x10 4)
SENA,=0
2.279 0.015
3.962 0.005
13.869 0,450
E 11,2=10-3
3.526 + 0.000
3.954 + 0.021
5.016 + 4.120
Ê 2 10-5
2.352 + 0.028
3.951 0.005
1.915 = 0.127
SENA J 1-0-3
1.936 0.023
3.967 0.023
15.211 7.351
SENA
1.103 + 0.010
3.938 + 0.022
15.978 9.905
SENA
1.009 + 0.004
3.951 0.009
12.492 3.054
MLP
1.036 0.012
3.962 0.019
3.872 0.225
Table 4: AE and VAE-based evaluation for SENA-6 across seeds. Models are sorted by sparsity.
AE-based
VAE-based
Method
Test MSE (x10-
Test MSE (x10-
Test DRL (x10-4)
SENA-ôx_0
2.252 0.011
3.951 0.012
0.007 = 0.009
11,x=10-
3.537 0.098
3.944 0.007
5.742 = 3.772
4 1-1==105
2.581 + 0.011
3.970 0.006
5.351 + 6.089
SENA-8
1.552 + 0.012
3.973 + 0.004
0.001 + 0.001
SENA- 11AA--x--000 0.01
1.096 + 0.020
3.969 0.027
0.024 + 0.005
SENA-8 0.1
1.012 +.000
3.966 0.010
0.326 0.111
ML.P
1.350 0.029
3.954 0.029
3.816 + 1.331
20

---

Published as a conference paper at ICLR 2025
V BENCHMARKING ON THE WESSELS DATASET
We included a second large-scale Perturb-seq dataset based on CRISPR-cas13 which aims at efficiintyly largeting multiple genes for combinatorial perturbations Wessels al. (2022). This technique, termed CaRPool-seq, encodes multiple perturbations on a cleavable CRISPR array that is associated with a detectable barcode sequence. CaRPool-seq was applied to THP1 cells, an acute myeloid leukemia (AML) model system, to perform combinatorial perturbations of myeloid differentiation regulators and identify their impact on AMI. differentiation phenotypes.
The perturbations include 28 single perturbations, 26 regulator genes and two negative control genes, as well as 158 double-gene perturbations. We performed standard preprocessing for single cell data (filtering, normalization, and log traansffroommion), yielding to total of 424 unperturbed cells (controls), 3036 cells under the 28 single-gene perturbations, and 17592 cells with double-gene perturbations. Pseudo-bulk expression profiles are then obtained by adding gene expression of cells sharing the same perturbation. The resulting profiles are then visualized as UMAP landscape. The same procedure performed for the Norman2019 dataset (Fig. 11 A and B). Marker genes then obtained using the Wilcoxon test to contrast the cells from each perturbation with the remaining cells implemented in Seurat's F Findll 1Ma rke rs, requiring adjusted p-value 0.001 and log-fold change> > 2.
The Wessels2023 study focused on perturbing myeloid differentiation regulators Wessels et al. (2022). This resulted in all perturbations having similar effects at the transcriptomics levels, as shown in Fig. 11. While in the Norman2019 datasets cells affected by different perturbations tend cluster separately (panel A), most of the interventions in Wessels2023 are grouped together (panel B), indicating similar profiles. Moreover, the number of genes that are differentially expressed following perturbation is generally lower in the Wessels2023 study than in Norman2019 109 (Fig. 11-C), indicating that overall the perturbations the Wessels2023 data had a more limited effect.
The peculiarities of the Wessels2023 data lead to notable results: both the SENA-discrepan VAE and the discrepancy- VAE consistently assign single-gene perturbations to single latent factor (figure not shown). This can be interpreted as the models recognizing that all single-gene perturbations have similar effects.
VAA
In terms of predictive capabilities (Table 5), we observe that MMD performances the double perturbations are comparable only when considering higher values of the À parameter. This indicates that for this dataset, introducing interpretability in addition to representational capabilities is more difficult, possibly due to the peculiarities of the performed experiments.
21

---

Published as a conference paper at ICLR 2025
A
Norman2019
B
2ess
EP300
VPRBP
FL1
MYB
SUPT
KDM1A
INTS1
HDAC3
MED
GFI
SMARCD
Parstelooo
Norvadtte
Figure 11: Comparison of Norman and Wessels datasets. A. UMAP representation of single-gene and combinatorial perturbations captured in the Norman2019 dataset. Each point represents the pseudo-bulk expression profile of a genetic perturbation. B. UMAP representation of the Wessels2023 dataset similar A. Perturbations which share one gene and form clear clusters are shown the same color to enhance clarity, while the remaining perturbations are shown individually.C. Cumulative density distribution of perturbations by their respective number of marker genes (differentially expressed genes with Pady 00,001, average LFC > 2) per perturbation in the Norman2019 and Wessels2023 dataset.
Table 5: Performance comparison between SENA-dixcrepancy- VAE and discrepancy- VAE the Wessel2023 dataset across different lambda values and latent factors for double perturbation samples. Note that KLD and L1 losses not dependent the samples, but computed after the training process is finished.
Latent Dimension
Encoder
Metric
50
28
14
MMD4
0.1810
0.0001998
0.1768
0.0001898
0.1886 + 0.0001225 225
MSE
0.0767
+ DOUR0051
0.0840 + 0.0000331
0.0934 + 0.0000015
Original MLP
KLD
0.0137
0.000011
0.0142
0.0150 + 0,0000002
L14
0.0054
0.000000]
0.0028
0.0000001
0,0016 1 0.00000004
MMDI 0.2285
0.0011976
0.1904
0.0001348
0.1753
0.0000720
SENA- -
MSEL
0.0934
0.0000111
0.0849
0.0000108
0.0782
0.000000S
KLD
0.0130
1 0.0000022
0.0135
1 0.0000025
0.0107
1. 0.0000003
L1J
0.0013
0.0021
0.0000018
0.0052 1
MMDI 0.3637 0.002126
0.4267 + 0.0027825
0,4049 +
SENA-82-0.
MSEL
0.1002 0.0000012
0.1081
0.0000056
0.1086
0.000000S
KLD
0.0040 = 0,00000]]
0.0049
0.0000001
0.0066 + 0.0000015
L11
0.0090 + 0.0000211
0.0026 + 0.0000001
0.0013 + 0.0000001
MMD4
0.5675 0.0091299
0.6156
0.0041942
0.5725 + 0.0096761
MSF.
0.1 1088 + DO0RD037
0.1068 + 0.0000076
0.1131 + 0.0000011
SENA-6x
KIDD..
0.0019 +
0.0025
0.0027 +
1.14
0.0120 DO00R0014
0.0040
0.000001D
0.0014
0X0000011
22

---

Published as a conference paper at ICLR 2025
VI EVALUATING GEARS ON PREDICTING UNSEEN DOUBLE PERTURBATIONS
We trained GEARS, following the authors' recommendations, for 20 epochs on all single-gene perturbations from the Norman2019 dataset, and predicted the same double-perturbatiion used to evaluate MeNA--iicrepaacc--A.. and the original Jacrepancy-VAA. We repeated this experiment with different initializations.
We then computed the MMD and MSE between the predicted and the true double cuuble-perturbations gene expression profiles by forwarding randomly-selecter
unperturbed cells (control) and averaged the metrics across perturbations. We show the MMD in the main benchmarking (Table 1) and in this section's is Table compare against MSE. When analyzing these results, we found that the MSE on double perturbations exhibited scores 10 times larger than the one reported during training for validation folds (0.00368 0.000363 across different latent dimension sizes for one seed), which suggest potential lack of generalization. Moreover, MMD also showed significantly larger scores compared to SENA or the standard discrepancy- VAE. This may be justified by the fact that the MMD not used as loss function in GEARS, hence true and predicted distributions may differ significantly when MSE is large enough. Finally, in order to have a baseline within GEARS, we computed the MSE and MMD between the true expression of double perturbations and true expression of unperturbed cells, yielding an average value of MMD = 14.39 and MSE 0.085.
Overall, these results highlight that i) the hidden dimension size of GEARS is largely dependent of the prediction performance. ii) GEARS outperforms the defined baseline in terms of MSE for some range of latent dimensions. iii) GEARS fails to capture the underlying distribution of double perturbations, possibly due to the lack of distribution-dittan distancc loss function during training.
Table 6: MMD and MSE of GEARS on Norman2019's dataset for double perturbation prediction across several latent dimensions
Metric
Latent Dimension
105
70
35
10
MMDJ
14.9420
12.6036
0,345302
13,2590
0.589124
12.9774
0.751333
15.3099 1.
MSFL 0.0986 + 0009157
0.0446 + 0.18200
0.0562 + 0.01
0.0533 + 0019782
0.1131 + 0.00568
23

---

Published as a conference paper at ICLR 2025
VII
ADDITIONAL MATERIALS ON NORMAN2019 DATASET
A
Non-mesked gradiants
106
Frequency
109
Masked
105
10
101
10
B
-2.0
-1.5
-1.0
-0.5
0.0
0.5
10
105
10
105
101
10
C
-0.4
-0.2
0.0
0.2
0.4
0.6
0.001
1.000
105
106
103
105
101
10
D
-0.4
-0.2
0.0
0.2
0.4
0.01
1.00
105
108
103
105
10
10
-0.4
-0.2
0.0
0.2
0.4
0.1
1.0
Figure 12: SENA layer analysis. Masked and mon-masked gradients for SENA SENA-discrepancy- VAE at the output of the SENA layer, for values of (A), 10 (B), 0.01 (C) and 0.1 (D), Barplot showing matrix M values is depicted next to each histogram.
24

---

Published as a conference paper at ICLR 2025
LHX1
TBX3
SPI1
FOXF1
COL1A1
EGR1
FOXL2
10"
10
10"2
10
10"
10
10"
10"4
False True
Falee
Falee
True
True
Falee
Fale
True
Affected CEBPA
CEBPB
CITED1
JUN
HOXA13
FOXA1
UBASH3B
10
10
10"2
10
10
10"2
10 :
10
10
10
10"
10
False
True
Falee
True
Faloe
True
Faloe
True
False
True
Falee
True
Falos
True
SGK1
TBX2
TGFBR2
TMSB4X
SLC4A1
PTPN12
POU3F2
10
10"22
10
10-2
10
10
10
10
10'9
1
10
10-2
10"4
10"
Faise
Thue
False
True
Faiso
Thue
False
Thuo
Fatse
True
Faise
Trua
False
Thuo
HNF4A
CDKN1C
HOXC13
LYL1
MEIS1
OSR2
ARRDC3
10
10
10
10
10-2
10-2
10"
PRDM1
ETS2
COL2A1
CNN1
CLDN6
CKS1B
CDKN1A
10
10
10-2
10
10 3
1
10
10
10
10
10
False
True
Falee
True
Falee
True
False True
Falze
True
Falee
True
False True
BCL2L11
AHR
10
10"
Figure 13: Differential Activation Scores analysis. DA analysis for targeted and non-affected BPs along the 37 single-gene perturbations present in the input gene expression matrix.
25

---

Published as a conference paper at ICLR 2025
Figure 14:
Mapping between latent factors and perturbations genes. Mapping dis-
tribution of knocked out genes and latent factors, from SENNA-disreppacc-WWe -AA (A = 0) and Aiscrepancy-VAF (MLP encoder) for several values of latent dimensions (105, 70, 35, 10, 5). We generated this mapping from the interventional encoder according to h = Softmax (Linear (LeakyReLU (Linear(c))) temp), where is the one-hot encoding vector for each perturbation and temp was set to 100 as recommended for inference in the original manuscript.
Table 7: Causal graph (Fig. 2) details. Each row lists latent factor, the number of targeted perturbations, and associated biological processes within
Latent Factor
Targeting Perturbations
Biological Processes
41
41
57
65
6
10
18
14
53
25
10
69
1
1
12
9
10
Table Top 6 (knockout, gene set) pairs according to DA score. Fig. details). DA scores for selected genes and associated GO terms.
Gene
GO Term
DA Score
COLIAI
GO:0038065
0.92
TBX3
GO:0006833
0.88
GLBIL2
GO:0070141
0.71
TP73
GO:0006833
0.69
CEBPA
GO:0010714
0.65
CEBPA
GO:0006595
0.60
26

---

Published as a conference paper at ICLR 2025
Table 9: Mapping between BP and latent factors.
Latent Factor
GO ID
GO Term
2
GO:2000352
negative regulation of endothelial cell apoptotic process
GO:0010944
negative regulation of transcription by competitive promoter binding regulation of ERAD pathway
2
GO:1904292
2
G0:0060216
definitive hemopoiesis
GO:0048557
embryonic digestive tract morphogenesis
2
GO:0071549
cellular response to dexamethasone stimulus
GO:0055023
positive regulation of cardiac muscle tissue growth
2
GO:0071243
cellular response to arsenic-containing substance
GO:0006067
ethanol metabolic process
2
GO:0099188
postsynaptic cytoskeleton organization water transport
GO:0006833
2
GO:0045723
positive regulation of fatty acid biosynthetic process regulation unsaturated fatty acid biosynthetic process ether metabolic process response to UV-A
GO:2001279
2
GO:0018904
12
G0:0070141
12
GO:0060512
prostate gland morphogenesis positive regulation of viral life cycle phasic smooth muscle contraction
12
GO:1903902
12
GO:0014821
12
GO:0006359
regulation of transcription by RNA polymerase III regulation of granulocyte differentiation
12
GO:0030852
12
GO:0044849
estrous cycle
12
GO:0045737
positive regulation of cyclin-dependent protein serine/threonine kinase activity
12
GO:1903319
positive regulation protein maturation
12
GO:2000010
positive regulation of protein localization cell surface hydrogen peroxide biosynthetic process response to lithium ion positive regulation of cellular extravasation
15
G0:0050665
41
GO:0010226
GO:0002693
GO:0045654
positive regulation megakaryocyte differentiation uterus development
GO:0060065
GO:1902042
negative regulation of extrinsic apoptotic signaling pathway via death domain receptors
GO:0000305
response to oxygen radical
41
G0:0001915
negative regulation of T cell mediated cytotoxicity
GO:0002281
macrophage activation involved in immune response defense response tumor cell
GO:0002357
GO:0006595
polyamine metabolic process
GO:0006921
cellular component disassembly involved in execution phase of apoptosis
GO:0006957
complement activation, alternative pathway
GO:0010714
positive regulation of collagen metabolic process
GO:0010829
negative regulation of glucose transmembrane transport
GO:0014912
negative regulation of smooth muscle cell migration cytolysis
GO:0019835
GO:0030449
regulation complement activation
GO:0032703
negative regulation of interleukin-2 production transforming growth factor betal production positive regulation of mast cell activation
GO:0032905
G0:0033005
GO:0039532
negative regulation cytoplasmic pattern recognition receptor signaling pathway
GO:0044342
type B pancreatic cell proliferation adhesion of symbiont to host transcytosis
41
GO:0044406
GO:0045056
GO:0045663
positive regulation of myoblast differentiation negative regulation of calcium-meding signaling actin crosslink formation
GO:0050849
GO:0051764
41
G0:0061450
trophoblast cell migration
GO:0070102
interleukin Beceleukki---medittee 6-mediated signaling pathway
41
GO:0070486
leykocyte aggregation
GO:0071391
cetlular response to estrogen stimulus

---

Published as a conference paper at ICLR 2025
Table 10: Mapping between BP and latent factors (Cont.).
Latent Factor
GO ID
GO Term
41
GO:0097202
activation of cysteine-type endopeptidase activity
41
GO:0099515
actin filament-based transport
GO:0140374
antiviral innate immune response
41
G0:1900120
regulation of receptor binding
GO:1900745
positive regulation of p38MAPK cascade
41
G0:1902307
positive regulation of sodium ion transmembrane transport regulation dendritic cell chemotaxis
GO:2000508
41
GO:0020027
hemoglobin metabolic process innate immune response mucosa tooth mineralization
GO:0002227
41
GO:0034505
GO:0038065
collagen-activant signaling pathway
GO:0001958
endochondral ossification
GO:0034389
lipid droplet organization
G0:0010715
regulation of extracellular matrix disassembly
G0:0051131
chaperone-mediate protein complex assembly innervation
GO:0060384
G0:0061684
chaperone-medine mediated autophagy
GO:0055091
phospholipid homeostasis
GO:0060742
epithelial cell differentiation involved in prostate gland development
GO:0097066
response to thyroid hormone
GO:0042976
activation of Janus kinase activity cell proliferation involved in heart morphogenesis
GO:0061323
41
GO:0002713
negative regulation of B cell mediated immunity sequestering of metal ion regulation programmed necrotic cell death regulation of cAMP -eMP-dependent protein kinase activity (eeduum-indeeendeen cell-cell adhesion via plasma membrane cell- adhesion molecules
GO:0051238
GO:0062098
41
GO:2000479
53
GO:0016338
53
GO:0009713
catechol-comtaining compound biosynthetic process bleb assembly
53
GO:0032060
53
GO:0036035
osteoclast development
53
G0:0007190
activation of adenylate cyclase activity
53
GO:0009065
glutamine family amino acid catabolic process endothelial cell chemotaxis
53
GO:0035767
53
GO:2000678
negative regulation of transcription regulatory region DNA binding positive regulation of macrophage activation
53
GO:0043032
53
GO:0060749
mammary gland alveolus development
65
G0:0021516
dorsal spinal cord development protein nitrosylation
65
GO:0017014
65
GO:0010842
retina layer formation
65
GO:0030540
female genitalia development
65
GO:0042759
long-chain fatty acid biosynthetic process
65
G0:2000696
regulation of epithelial cell differentiation involved in kidney develop- ment
65
G0:0015701
bicarbonate transport
65
GO:0034638
phosphatidylcholle catabolic process lymph node development
65
GO:0048535
65
GO:0050995
negative regulation of lipid catabolic process endothelial cell morphogenesis
69
GO:0001886
28


Reviews:
[{"review_id": "RCSrbP35wp", "reviewer": "Reviewer_Gg4s", "summary": "The paper addresses the challenge of learning interpretable causal representations for Perturb-seq data (gene expression in cells). The primary contribution is the novel introduction of masking to incorporate biological process (BP) knowledge into an existing method for causal representation learning (discrepency-VAE), which is named SENA-discrepancy-VAE. The masking ensures that latent factors can be interpreted as linear combinations of the activity of BPs. Since this modification is compatible with the discrepancy-VAE, the original model's theoretical guarantees for causal representation learning remain.\n\nThe method and ablated variants are evaluated on a Perturb-seq dataset collected from one particular cell line and is set up to minimize the overlap between the BPs. The results demonstrate that SENA performs similarly to discrepency-VAE in terms of reconstruction yet results in sparser and more interpretable results. Furthermore, by studying the contrast between inferred activity levels on perturbed and control samples the authors show that the latent factors can be associated with BPs and are therefore interpretable.", "questions": "- How sensitive is the model to the quality and completeness of the pathway knowledge used? Have you tested with different pathway databases or subsets of pathways?\n- How does the computational complexity scale with the number of biological processes? Is there a practical limit to how many processes can be incorporated?\n- Have you explored whether the causal relationships discovered by the model align with known biological pathway interactions beyond the examples provided?\n- Can you confirm that all genes involved in the double perturbations were also present in your single-perturbation training data?\n- Has $N$ and $\\tau$ been mixed up in the Hits@N metric?\n- How does table 2 show that “both models tend to assign most interventions to a small number of latent factors”?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "- Clear technical contribution that bridges causal representation learning with biological interpretability while maintaining theoretical guarantees\n- The paper contributes to causal representation learning for Perturb-seq data by introducing biological interpretability through pathway information, while maintaining the theoretical guarantees of discrepency-VAE.\n- Well written and clear presentation of the method and results.\n- Thorough ablation studies\n- Demonstrates interpretability of latent factors with concrete biological examples.", "weaknesses": "- Experimental validation is limited to one dataset and no baselines other than their own ablations and discrepancy-VAE. The paper would benefit from comparisons to at least one of the other listed related works.\n- No comparison with simpler approaches like post-hoc interpretation of standard discrepancy-VAE latent factors.\n- While the link between latent factors and BPs is investigated, the quality of the discovered causal graph is not.\n- Given that the latent factors group a large number of BPs into a small number of latent factors there should be a deeper investigation of the biological plausibility and practicality of this result beyond the contrasting activations.\n- Readability of several figure texts should be improved.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "2XwCciqoZv", "reviewer": "Reviewer_YiHQ", "summary": "This work proposes to extend the discrepancyVAE interventional causal representation learning framework to biological processes applications. Specifically, the authors propose to embed prior knowledge about biological processes (BPs) through a framework called SENA-discrepancyVAE, which recovers latent factors that are a linear combination of a set of biological processes (pathways). The main idea presented in this work is to design a more flexible encoder class (SENA-\\delta) specific to mapping biological pathways to latent causal factors for interpretability. Empirical results show that the framework is shown to improve performance in predicting the effect of unseen perturbations.", "questions": "What is the intuition behind the $\\lambda$ hyperparameter to tune small influences of a gene on a biological process? Should this be a constant value throughout the mask matrix or would it be better to learn this influence via some type of attention weights?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- The paper is well written with strong motivations behind using CRL techniques for biological applications.\n- The metrics proposed (differential activation, Hits@N) seem to be robust indicators of perturbation effects on BPs and downstream effects. I believe these evaluation metrics are one of the key interesting contributions of this work.\n- The empirical evaluation is exhaustive and illustrates some interesting observations, especially the representational capacity of the VAE-based SENA method compared to the traditional discrepancyVAE.\n- The interpretability analysis of the reparameterization layer is interesting and reveals which genes were affected the most upon perturbations. I do believe that exploring real-world applications of CRL is a very important direction.", "weaknesses": "- Although the application in gene regulatory networks is quite interesting, this work seems to be more of an evaluation study of the discrepancy-VAE framework proposed by Zhang et al. I do not see much of an added contribution beyond the original paper besides highlighting the application.\n- The difference in performance between the SENA variant and the original discrepancyVAE seems to be quite marginal in terms of representation in the double-perturbation scenario. For instance, in Table 2, the KL-divergence for double-perturbation prediction is only marginally better than the original MLP-based discrepancyVAE.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "CnUimR2dcr", "reviewer": "Reviewer_ez4D", "summary": "This paper presents SENA-discrepancy-VAE, an extension of the discrepancy-VAE framework that incorporates biological pathway knowledge to produce interpretable causal latent factors. The authors modify the encoder architecture to map gene expression through biological processes (BPs) while maintaining the theoretical guarantees of the original model. The approach achieves comparable predictive performance to the original discrepancy-VAE while providing biologically meaningful latent representations and interpretability.", "questions": "Can you provide examples of BPs in the appendix? \n\nRegarding DAR evaluation: What happens if unnaffected pathways have very low action? Also, more general, how would you deal with impalanced pathways, which might lead in measuring large noise levels? \n\nIn table 2, for SENA λ=0.1, latent dim 105, the variance compared to original MLP and λ=0 is significantly lower (0.000081 vs 0.001087). Just double checking if this is corrent.\n\nL237: During filtering you end up with a (biased) set of BPs. How much do you think this can influence interpretability? Is there a risk of removing useful BPs?\n\nSuggestions\nL100: the word faithfully here gets confused with causal faithfulness. Please consider an alternative adverb if possible.\nL105: target instead of targets (remove final s)", "limitations": "", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "* Clarity: The paper is well written and easy to follow.\n\n* Novel Technical Contribution: The paper successfully extends causal representation learning to incorporate domain knowledge while preserving theoretical guarantees. The SENA-δ encoder architecture is a clever solution to balance interpretability and performance.\n\n* Practical Impact: The work addresses a significant gap in current causal representation learning methods for biological data, where interpretability is crucial for scientific insights. The ability to map latent factors to biological processes makes the model more useful for domain experts.", "weaknesses": "* Limited Biological Validation: While the authors show statistical associations between perturbations and biological processes, there could be more validation using external biological knowledge or experimental validation of the discovered causal relationships.\n\n* Hyperparameter Sensitivity: The model introduces an additional hyperparameter λ that significantly impacts performance. While ablation studies are provided, more guidance on selecting this parameter would be valuable (this is important given that there's some large impact on the performance of the method)\n\n* Restricted Evaluation: The empirical evaluation is limited to a single dataset (Norman et al., 2019). Additional validation on different types of biological data would strengthen the claims of generalizability.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "PJp2j6lXB0", "reviewer": "Reviewer_ox26", "summary": "This paper introduces SENA-discrepancy-VAE, a novel model in causal representation learning (CRL) designed to make biological data analysis—especially from Perturb-seq experiments—more interpretable. A key innovation of this work is how it integrates biological processes (BPs) as prior knowledge, directly linking the model’s latent factors to known biological pathways. This approach fills an important gap in existing CRL methods, which often struggle with interpretability since they don't directly associate learned representations with actual biological mechanisms, making them less useful for real research applications.\n\nSENA-discrepancy-VAE builds on the standard discrepancy-VAE by introducing a pathway-based masking strategy within a new encoder, SENA-δ. This encoder uses a two-layer masked MLP where the first layer maps gene expression values to BP activity levels, with a tunable parameter that adjusts the influence of genes outside predefined pathways, giving the model flexibility in gene-pathway associations. The second layer models latent factors as combinations of these BP activities, which is a more realistic approach since biological interventions often impact multiple pathways. This setup stays true to the CRL assumption that each intervention targets a single latent factor but does so in a way that aligns with biological realities.\n\nThe authors evaluate SENA-discrepancy-VAE on a Perturb-seq dataset of leukemia lymphoblast cells. They show that the model performs as well as the original discrepancy-VAE on unseen perturbation combinations while providing greater interpretability by identifying specific BPs associated with each latent factor. This interpretability is validated through pathway-specific analysis, demonstrating the model’s ability to reveal biologically meaningful patterns in response to genetic interventions.", "questions": "1) Could the authors clarify whether multi-label or probabilistic pathway activity labels were considered as an alternative to binary labels? Binary labels may oversimplify gene activity levels, especially when certain pathways exhibit gradations rather than discrete on/off states.\n\n2) How would the model perform if pathway relevance were dynamically adjusted based on specific task contexts or experimental conditions? Pathway importance often varies, and adapting pathway relevance could improve model flexibility.  A task-specific analysis to explore whether dynamically adjusting pathway selection enhances generalizability across datasets and biological contexts could provide insights into the model’s robustness.\n\n3) In the causal graph (Figure 3), how sensitive are the edge weights to the choice of λ and latent dimension? Please provide a sensitivity analysis.\n\n4) The mask matrix M (equation 2) assumes binary gene-pathway relationships. Have you considered using weighted relationships based on pathway databases' confidence scores?\n\n5) How does D_KL evolve during training for different λ values? Training curves would help understand if this is an optimization or regularization effect. Does this phenomenon persist if you randomize the pathway annotations while maintaining the same sparsity structure? This would determine if the benefit comes from biological knowledge or just sparsity.\n\n6) Have the authors empirically validated that the expectation in Equation (10) aligns with the observed experimental outcomes? Demonstrating this match would reinforce the theoretical assumptions and provide additional confidence in the model's causal interpretability.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1) The proposed model’s integration of biological pathway data as a prior in the causal representation learning (CRL) framework is interesting and practical, offering a biologically grounded solution to gene expression analysis. Through the direct alignment of the latent factors with biological pathways, SENA-discrepancy-VAE addresses the common limitation in CRL models of producing uninterpretable latent factors. \n\n2) The paper presents thorough experiments across multiple perturbation types, including single and double-gene knockouts, demonstrating the model’s robustness. Its generalization to unseen perturbations is compelling, and the ablation studies, which explore interpretability-reconstruction trade-offs, further validate the model’s design.\n\n3) The authors have done a good job communicating the importance of embedding biological processes into latent spaces, with visuals that illustrate BP-specific activity levels influencing latent factors. The use of causal graphs and the differential activation (DA) metric enhances transparency, allowing readers to trace latent factors back to biological functions. \n\n4) The model’s ability to provide interpretable predictions on cellular responses can aid in experimental design and offer insights into the potential effects of genetic or drug interventions. This approach addresses a pressing need in biomedicine for interpretable causal representation learning (CRL) models that can shed light on the intricate causal relationships underlying gene function and cellular processes.", "weaknesses": "1) While the model effectively identifies single-point perturbations, it does not accommodate multi-step perturbations or capture the progression of cellular responses over time. In biological experiments, the cellular response often evolves in phases, with gene activity showing distinct transitions that are essential for understanding the effect of interventions.\n\n2) The model’s validation on a single dataset (K562 cell line data) restricts insights into its generalizability across different cell types or conditions. Testing on additional datasets, such as those from other cell lines or cellular environments, would offer a stronger assessment of robustness and applicability across a wider range of biological data.\n\n3) The model assumes static pathway relevance across all tasks, which may limit its adaptability in varied biological contexts where pathway importance changes with cell type or condition.\n\n4) The paper assumes that each intervention corresponds to a single latent factor, limiting the model's ability to capture complex interactions where multiple latent factors might be influenced by a single intervention. This simplification restricts the model’s interpretability in representing overlapping or interacting biological processes, which are common in gene expression dynamics.\n\n5) The model lacks a detailed exploration of how varying the number of latent dimensions impacts the interpretability and causal mapping of latent factors. Larger or smaller dimensions can influence the granularity of the factors and thus affect the biological insights the model can provide.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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:

Decision: Accept (Poster)