Automated Discovery of Pairwise Interactions from Unstructured Data

## Abstract

Abstract Pairwise interactions between perturbations to a system can provide evidence for the causal dependencies of the underlying underlying mechanisms of a system. When observations are low dimensional, hand crafted measurements, detecting interactions amounts to simple statistical tests, but it is not obvious how to detect interactions between perturbations affecting latent variables. We derive two interaction tests that are based on pairwise interventions, and show how these tests can be integrated into an active learning pipeline to efficiently discover pairwise interactions between perturbations.
We illustrate the value of these tests in the context of biology, where
pairwise perturbation experiments are frequently used to reveal interactions
that are not observable from any single perturbation.
Our tests can be run on unstructured data, such as
the pixels in an image, which enables a more general notion of interaction than
typical cell viability experiments, and can be run on cheaper experimental assays.
We validate on several synthetic and real biological experiments that our tests are able to identify interacting pairs effectively.
We evaluate our approach on a real biological experiment where we knocked out
50 pairs of genes and measured the effect with microscopy images. We show that
we are able to recover significantly more known biological interactions than
random search and standard active learning baselines.

### 1 Introduction

Figure: Figure 1: The traditional experimental design loop involves a number of human expert decisions: the expert need to *measure* a carefully chosen feature of the experimental outcome (e.g. did a cell survive a perturbation?), formulate a prediction for the behaviour of this feature under the null hypothesis to test for interactions, and select interacting pairs for the combinatorial space of possible perturbations. Our approach enables automated selection of perturbants and testing for interactions directly from raw signal data (e.g. images of cells under a microscope).
Refer to caption: x1.png

Across the sciences, measurement of pairwise interaction between perturbations often reveals the existence of underlying mechanisms that single perturbations cannot. For example, quantum entanglement experiments support the counterintuitive predictions of quantum physics by demonstrating that entangled particles’ spins which we would expect (under classical laws) to be independent, are in fact perfect anti-correlated once we make a measurement.
In economics, people may be more (or less) willing to pay for goods when presented in a bundle than they are to pay for each good in isolation, which reveals complements (or substitutes) in the underlying consumer preferences.
And in biology—the application area on which we focus—the concept of *synthetic lethality* (Nijman, ) occurs when the simultaneous perturbation of two genes results in cell death while individual perturbations do not. Synthetic lethality reveals that the genes’ associated proteins play redundant roles in the underlying cellular mechanisms.

In order to demonstrate a pairwise interaction a scientist (the human expert) has to carefully design three steps of the experiment:

- 1.
Measurement: the human expert to decides on what specific properties of the system to measure in order to reveal the interaction. For example, measuring a particle’s spin may reveal entanglement, while the mass of a particle will not; similarly measuring cell viability may reveal synthetic lethality, while measuring the cell’s colour will not.
- 2.
Hypothesis testing: an interaction is, by definition, a deviation from the effect that we expect under some null hypothesis which assumes each perturbation acts independently. The expert needs to specify the expected outcome assuming independence and compare this prediction to the actual outcome.
- 3.
Selection: there are typically many variables that one could perturb (e.g. there are roughly 20 000 genes in the human genome), but only a small fraction of them will exhibit surprising interactions. The expert selects pairs of perturbations from this space of all possible pairs.

Measuring interactions is further complicated by the fact that measurement and hypothesis testing are interdependent. We can see this in synthetic lethality example, where if we choose to measure the fraction of cells that survived the respective perturbations, then in order to test independence we would need to test whether $P(\text{survive}_{a}\cap\text{survive}_{b})=P(\text{survive}_{a})P(\text{
survive}_{b})$.
If we had instead measured the fraction of cells that died, we would have needed to test whether $P(\text{die}_{a}\cup\text{die}_{b})=P(\text{die}_{a})+P(\text{die}_{b})-P(
\text{die}_{a})P(\text{die}_{b})$. The state of the cells in the respective petri dishes is the same in both cases, but correctly testing (in)dependence depends on how the expert chose to measure that state. These three steps require significant expertise and knowledge to choose and conduct the measurements, which makes scientific discovery hard to automate and scale.

Modern high throughput screening platforms (Dove, ; Baillargeon et al., ; Blay et al., ; Bock et al., ; Fabio et al., ; Morris et al., ; Messner et al., ) allow us to run large scale perturbation experiments that collect information-rich unstructured measurements. For example, cell painting assays (Bra, ; Fay et al., ; Chandrasekaran et al., ), provide microscopy images of cells which capture the same cells that the human expert would have measured, but without pre-committing to measuring a particular property of the cell, such as whether it is alive or dead.
This lets us *measure* all the potential properties of interest at once—albeit as unstructured signal data, such as raw pixel images or sensor measurements, rather than preprocessed properties—but it is not clear how to use these measurements for interaction *testing* and *selection* in order to efficiently discover pairwise interactions.
This suggests the primary question that this paper seeks to address: if we dispense with prespecified structured measurements, can we automate the discovery of pairwise interactions by designing testing and selection procedures that efficiently detect interactions from unstructured data?

To test for interactions,
we show that pairwise perturbations are *separable* if the pairwise experiments provide no additional information beyond what is already known from the single perturbations.
This separability property can be tested by comparing density ratios of the single and double perturbations to the base / control distribution.
We then develop an analogous test for the case where pairs of perturbations affect disjoint subsets of the outcome space. For example, if two perturbations affect different organelles in the cell, each of which affect disjoint pixels in our images.
This strong form of non-interaction is of interest because it gives sufficient conditions for when we can compose summary statistics of the form $\mathbb{E}[h(x)]$ from single perturbations, to predict the corresponding statistics from double perturbations. In particular, this lets us compose embeddings of single perturbations to predict the embeddings of double perturbations in a manner analogous to the classic word-vector analogies Mikolov et al. ().

With this notion of an interaction, we address the second task of selection via efficient experimental design. We can search the space of pairwise experiments by selecting pairs of perturbations that are likely to result in large test statistics.
In doing so, we reduce the problem of finding interacting pairs of perturbations into an active matrix completion problem. By defining this notion of interaction, we avoid the need to characterize uncertainty over the pixel-level outcomes (or some embedding thereof). Instead we directly model the posterior over an unknown reward matrix defined by the test statistics and at every round select actions striking a balance between exploration and exploitation via information directed sampling (IDS; Russo and Van Roy, ; Xu et al., ).

We evaluated our approach on both synthetic and real biological experiments. On the synthetic tasks, we found that our tests were able to detect both forms of dependence, validating our theory. On a benchmark consisting of all pairs of gene knockouts for 50 genes in HUVEC cells, we found that our approach using IDS discovers pairs of genes with higher interaction scores significantly faster than baselines, resulting in $10\%$ more known biological interactions being discovered and potentially far more novel interactions. The interactions we detected were also complementary to those which would have been discovered using existing cosine similarity-based approaches, and as a result, the two approaches can be combined to get a more detailed estimate of the relationships between genes from perturbation experiments.

In summary, we develop a system for discovering pairwise interactions through the following contributions,

- •
We show that two perturbations are separable if and only if, double perturbations provide no additional information than that which was revealed in single perturbations.
- •
We show that two perturbations are disjoint if and only if their density compose additively. This result also implies sufficient conditions for embeddings of perturbations to be composed additively to predict pairwise perturbations.
- •
We show that the test scores for detecting interactions can be used to efficiently search for pairs of perturbations that interact using active matrix completion.

### 2 Related work

###### Causal representation learning

Our approach to the problem of detecting interactions builds on the modelling assumptions developed in nonlinear independent component analysis (Hyvärinen and Pajunen, ; Hyvarinen and Morioka, ; Hyvarinen et al., ) and causal representation learning (Schölkopf et al., ), where we assume we observe some nonlinear mixing function, $g(\cdot)$, of the underlying latent variables.
The causal representation learning literature typically focuses on disentangling latent variables, whereas we look for testible implications without disentanglement. Our data generating process assumptions are most similar to the interventional setting (Ahuja et al., ; Squires et al., ; Buchholz et al., ; Varıcı et al., ).
If we could successfully disentangle latent variables, then our task would be straightforward, but disentanglement is challenging in practice because it is not possible to validate whether an algorithm has succeeded in disentangling latent variables without access to ground truth. Like us, Jiang and Aragam () also attempt to learn latent dependence properties (they characterize the whole latent DAG, not just marginal dependence) without disentanglement, but they assume a fixed bipartite graph of dependence between latents and observations, which does not apply to the pixel-level observations that we study.
Zhang et al. () develop causal representation learning techniques to disentangle latent variables and characterize conditions for extrapolation. Like us, they focus on biological applications, but they rely on stronger polynomial assumptions to achieve disentanglement. Their conditions for extrapolation are dual to our separability tests in that they argue that you can extrapolate when interventions affect *non-overlapping* latents, while we seek to discover when intervention affect *overlapping* latent.
Finally, our separability test is closely related to, and inspired by Wang et al. (), but whereas they assume separability of concepts to manipulate generative models, we aim to test for an analogous notion of independence in experimental data.

###### Representation learning for gene knockouts

Our experiments build on a number of recent works showing the effectiveness of embedding cells presented into a representation space Sypetkowski et al. (); Kraus et al. (); Xun et al. (). These works show that you can infer that genes code for proteins that form part of complexes by finding embeddings that are highly cosine similar. This works because knocking out proteins on the same pathway will induce the same morphological effect, but this approach is limited to effects that are revealed by single perturbations.
There have also been a number of papers that have attempted to learn disentangled representations of cells (e.g. Lotfollahi et al., ; Bereket and Karaletsos, ; Lopez et al., ), mostly from gene expression data.

###### Design of Gene Knockout Experiments

There has been considerable work in the recent years to develop algorithms for the design of gene knockout experiments. Typically the problem is studied in the context of discovering single gene-knockouts (Mehrjou et al., ) which result in a particular phenotypic effect of interest. A variety of methods including bandits (Pacchiano et al., ) and traditional experimental design (Lyle et al., ) approaches have been studied in this context. Further, approaches which aim to learn predictors for the effects of unseen gene knockouts typically operate on RNA-seq data (Huang et al., ). RNA-Seq data is far more structured than the image data we consider: the former is a high-dimensional vector of count data where each element corresponds to a particular gene’s expression, while the latter is just an unstructured collection of pixel intensities.

###### Experimental design.

The task of designing of experiments that acquire the most information about the system efficiently can be formalized as sequential Bayesian optimal experimental design (BOED; Ryan et al., ; Foster, ; Rainforth et al., ), where the goal to design experiments $x\in\mathcal{X}$ with outcomes $y\in\cal Y$ governed by a generative process $y\sim p(y\mid\gamma,x)$ with parameters $\gamma$. The experiments are performed sequentially $(x_{1},\dots,x_{T})$, with the objective of maximizing a measure of utility: the information gain (Lindley, ; Sebastiani and Wynn, ).
While this framework is elegant, it is challenging to apply
when $x$ is high-dimensional because information gain comparisons are very difficult to estimate in high dimensional data. We avoid this by defining the explicit task of detecting interactions.

### 3 Tests for pairwise interactions

In this section we develop tests for both *separable* and *disjoint* interventions. The
separability test, Section 3.1, allows us to test whether two
perturbations act on disjoint sets of latent variables. Disjointedness,
Section 3.2, implies compositional generalization of summary
statistics, allowing us to reduce the search space by predicting the
outcome of experiments without explicitly running them.

###### Setup

We have observations $X$ in an *observation space* ${\mathcal{X}}$ of unstructured measurements such as pixels in an image,
and a finite set of perturbations $\{T_{i}\in{\mathcal{T}}_{i}:i\in[n]\}$.
Although our theory accommodates generic perturbations, we restrict the discussion to binary perturbations for convenience,
i.e., ${\mathcal{T}}_{i}:=\{0,1\}$ for all $i\in[n]$, and $T_{i}=1$ means perturbation $i$ is applied.
For all $i,j\in[n]$, denote the perturbation indicator as follows:

$$ $\displaystyle\delta_{0}$ $\displaystyle=\left\{T_{[n]}=0\right\}$ $\displaystyle\delta_{i}$ $\displaystyle=\left\{T_{i}=1,T_{[n]\setminus\{i\}}=0\right\}$ $\displaystyle\delta_{ij}$ $\displaystyle=\left\{T_{i}=T_{j}=1,T_{[n]\setminus\{i,j\}}=0\right\}.$ (2) $$

In this section, we assume that for a pair of perturbations $T_{i},T_{j}$, we have access to experimental data
from four distributions: $p(x|\delta_{0}),p(x|\delta_{i}),p(x|\delta_{j})$, and $p(x|\delta_{ij})$; when we discuss the active learning procedures in Section 4, we will assume that we have access to all single perturbations distributions, $p(x|\delta_{i})$, and we will adaptively select the pairs, $i,j$, on which to collect samples from $p(x|\delta_{i},\delta_{j})$.
Finally, we assume the existence of a set of latent variables
$\{Z_{1},\dots,Z_{L}\}\subseteq{\mathcal{Z}}$ in some latent space, ${\mathcal{Z}}$, that capture
all relevant information about the perturbation. We state this precisely as,

###### Assumption 3.1 .

$X\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 4.0mu{
\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 4.0
mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 4.
0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}
\mkern 4.0mu{\scriptscriptstyle\perp}}}(T_{1},\dots,T_{n})|Z$, or equivalently,

$$ $\displaystyle p(x|T_{1},\dots,T_{n})=\int_{{\mathcal{Z}}}p(x|z)p(z|T_{1},\dots ,T_{n})\mathrm{d}z.$ (3) $$

Throughout, we assume all distributions have well-defined densities or probability mass functions with respect to some $\sigma$-finite base measure on their corresponding sample spaces.
We use the same symbol to denote a distribution and its density.
All proofs are deferred to Appendix A.

###### Assumption 3.1 .

#### 3.1 Separability testing

It would be trivial to verify whether two perturbations intervene on different latents
if the underlying latent variables were observed directly.
However, our observations are unstructured signal data, which we assume is generated by some deterministic mixing function.

###### Assumption 3.2 .

There exists a diffeomorphism(^1^11A diffeomorphism is a differentiable bijection with a differentiable inverse.) $g:{\mathcal{Z}}\to{\mathcal{X}}$ such that $X=g(Z)$.

###### Remark 3.3 .

This diffeomorphism assumption between ${\mathcal{Z}}$ and ${\mathcal{X}}$ can be relaxed;
our theory and methodology remain valid as long as the change of variable formula holds for the distributions of $Z$ and $X$.
For example, the distribution of $X$ may have support on a low-dimensional manifold of ${\mathcal{X}}$, with the latent space ${\mathcal{Z}}$ having a lower dimension than ${\mathcal{X}}$.
See Krantz and Parks (, Lemma 5.1.4) for the generalized change of variable formula for non-bijective transformations.

The key observation that we leverage is that when latent variables are independent, the change of variable formula implies that the density ratio of a
perturbed distribution to the original (control) distribution has a simple form that only involves the distribution of the intervened *latent* variable,

$$ $\displaystyle\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{ eq:densityratio}}{e}q:densityratio}\frac{p(x|\delta_{i})}{p(x|\delta_{0})}= \frac{p_{Z}(g^{-1}(x)|\delta_{i})\left|\text{det}(J(g^{-1}(x)))\right|}{p_{Z}( g^{-1}(x)|\delta_{0})\left|\text{det}(J(g^{-1}(x)))\right|}=\frac{p^{\dagger}_ {Z_{i}}(g^{-1}(x))}{p_{Z_{i}}(g^{-1}(x))}.$ (4) $$

Here $p^{\dagger}_{Z_{i}}$ denotes the perturbed distribution of the latent
variable $Z_{i}$ targeted by intervention $i$. The testable conclusion that we can
derive from this observation, is that when two variables are independent, we can
predict the density ratio of the double perturbation,
$\frac{p(x|\delta_{ij})}{p(x|\delta_{0})}$, as the product of the density ratios of the respective single perturbations,
$\frac{p(x|\delta_{i})}{p(x|\delta_{0})}\frac{p(x|\delta_{j})}{p(x|\delta_{0})}$.
To make this rigorous, we assume that there exists a causal factorization of the
latent distribution, where each latent variable is conditionally independent of
its non-descendents given its parents. We use $\operatorname{Pa}(Y)$ (resp. $\operatorname{Ch}(Y)$) to denote the parent (resp. children) of
a random variable $Y$.

###### Assumption 3.4 .

The latent distribution $p_{Z}(z)$ can be factorized into $L$ latent factors:

$$ $\displaystyle p_{Z}(z|T_{1},\dots,T_{n})=\prod_{l=1}^{L}p_{Z_{l}}(z_{l}| \operatorname{Pa}(z_{l})),$ (5) $$

where $\forall l\in[L]$, $\operatorname{Pa}(Z_{l})\subseteq T_{[n]}\cup Z_{[L]}$,
and $\forall i\in[n]$, $\operatorname{Ch}(T_{i})\subseteq Z_{[L]}$.
The random variable $Z$ here should be interpreted as the concatenation of all latent factors $\{Z_{1},\dots,Z_{L}\}$.

The observed effects of two non-interacting perturbations are attributed to distinct causal pathways,
rather than being confounded by a shared latent factor.

###### Definition 3.5 .

Two perturbations $\delta_{i},\delta_{j}$ are separable if $\mathrm{Ch}(T_{i})\cap\mathrm{Ch}(T_{j})=\emptyset$.

When two perturbations act separably, the resulting density ratios of the observation distribution will have the predictable interactions that we derived above.

###### Theorem 3.6 .

Suppose that Assumptions 3.4 and 3.2 hold, and that $\delta_{i},\delta_{j}$ are separable. Then,

$$ $\displaystyle\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{ eq:separabldensityratio}}{e}q:separabldensityratio}\frac{p(x|\delta_{i,j})}{p( x|\delta_{0})}=\frac{p(x|\delta_{i})}{p(x|\delta_{0})}\frac{p(x|\delta_{j})}{p (x|\delta_{0})}.$ (6) $$

This provides a testable implication of separability, but to test it, we need to derive a real-valued test statistics. By taking logs and rearranging terms, we can rewrite Eq. 6, as testing,

$$ $\displaystyle\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{ eqn:log_rel}}{e}qn:log_{r}el}\log p(x|\delta_{ij})+\log p(x|\delta_{0})-\log p (x|\delta_{i})-\log p(x|\delta_{j})\stackrel{{\scriptstyle?}}{{=}}0$ (7) $$

If we instead take expectations of Eq. 7 with respect to $p(x|\delta_{0})$, then this
amounts to testing if

$$ $\displaystyle\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{ eq:kltest}}{e}q:kltest}\mathrm{D_{KL}}(p_{0}||p_{i})+\mathrm{D_{KL}}(p_{0}||p_ {j})\stackrel{{\scriptstyle?}}{{=}}\mathrm{D_{KL}}(p_{0}||p_{i,j}),$ (8) $$

where $p_{i}=p(x|\delta_{i}),p_{0}=p(x|\delta_{0})$ and $\mathrm{D_{KL}}(\cdot||\cdot)$ denotes the Kullback–Leibler (KL) divergence.
Note that KL divergence measures the difference between two distributions; in our context, it quantifies the distribution shift resulting from interventions.
Eq. 8 reflects the intuition that the influence of two separable perturbations, measured by KL divergence, is additive.
We use the KL score, $\left|\mathrm{D_{KL}}(p||p_{i})+\mathrm{D_{KL}}(p||p_{j})-\mathrm{D_{KL}}(p||p
_{i,j})\right|$, to quantify the violation of the separability of $\delta_{i},\delta_{j}$,
where these KL divergences are estimated using samples.

In practice, we observe that the simple K-nearest neighbor-based (KNN) estimator of
the KL divergence (Wang et al., ) suffices for
low-dimensional problems. However, for high-dimensional data such as images, more dedicated estimation procedures are required (Belghazi et al., ; Song and Ermon, ; Ghimire et al., ).
We detail our estimation procedure in Appendix B.

###### Assumption 3.2 .

###### Remark 3.3 .

###### Assumption 3.4 .

###### Definition 3.5 .

###### Theorem 3.6 .

#### 3.2 Disjointedness testing

Our second testing procedure examines whether two perturbations operate on disjoint domains such that their effects are cumulative.
For example, if the morphological changes from knocking out two non-interacting genes can be separated into distinct visual features, such that their measures sum, then they are disjoint. We can define disjointedness of two perturbations formally as,

###### Definition 3.7 .

Two perturbations $\delta_{i},\delta_{j}$ are disjoint if

$$ $\displaystyle\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{ eq:disjoint}}{e}q:disjoint}p(x|\delta_{ij})-p(x|\delta_{0})=(p(x|\delta_{i})-p (x|\delta_{0}))+(p(x|\delta_{j})-p(x|\delta_{0})).$ (9) $$

Disjointedness is important, because it implies that we can predict pairwise
summary statistics of the distributions from individual perturbations. To see
this, let $h(\cdot)$ denote *any* feature map from the observation, $X$.
If we have disjoint perturbations $\delta_{i},\delta_{j}$, then,

$$ $\displaystyle\mathbb{E}[h(x)|\delta_{i,j}]-\mathbb{E}[h(x)|\delta_{0}]$ $\displaystyle=\int h(x)[p(x|\delta_{ij})-p(x|\delta_{0})]\mathrm{d}x$ (10) $\displaystyle=\int h(x)[(p(x|\delta_{i})-p(x|\delta_{0}))+(p(x|\delta_{j})-p(x |\delta_{0}))]\mathrm{d}x$ (11) $\displaystyle=\mathbb{E}[h(x)|\delta_{i}]-\mathbb{E}[h(x)|\delta_{0}]+\mathbb{ E}[h(x)|\delta_{j}]-\mathbb{E}[h(x)|\delta_{0}].\addcontentsline{lla}{section} {\numberline{\string\crtrefnumber{eq:embeddingtest}}{e}q:embeddingtest}$ (12) $$

This implies that we can define average centered embedding vectors, $\vec{h}_{i}:=\mathbb{E}[h(x)|\delta_{i}]-\mathbb{E}[h(x)|\delta_{0}]$ and $\vec{h}_{j}$,
and accurately predict $\vec{h}_{i,j}=\vec{h}_{i}+\vec{h}_{j}$ without running the experiments.
It is worth noting that there is growing evidence in the literature (Lotfollahi et al., ; Gaudelet et al., ) that shows that these relationships often hold in real biological experiments (and our experiments support this).
Disjointedness explains the sufficient conditions for this to hold. Section C.2 explains the choice of $h$ we tested.

We can test whether Eq. 9 holds, by testing the following null hypothesis:

$$ $\displaystyle\mathrm{H}_{0}:\frac{1}{2}p(x|\delta_{ij})+\frac{1}{2}p(x|\delta_ {0})=\frac{1}{2}p(x|\delta_{i})+\frac{1}{2}p(x|\delta_{j}).$ (13) $$

Given interventional data from $p(x|\delta_{0}),p(x|\delta_{i}),p(x|\delta_{j})$ and $p(x|\delta_{ij})$,
we can frame this as a standard two-sample test problem.
With balanced experiments, i.e., $p(\delta_{0})=p(\delta_{i})=p(\delta_{j})=p(\delta_{ij})$,
we can create samples from the mixture $\frac{1}{2}p(x|\delta_{ij})+\frac{1}{2}p(x|\delta_{0})$ by combining data from the controlled and doubly perturbed groups.
Similarly, we can obtain samples from$\frac{1}{2}p(x|\delta_{i})+\frac{1}{2}p(x|\delta_{j})$.
When experiments are heavily unbalanced, we can balance the datasets through downsampling or upsampling.

In practice, we employ the maximal mean discrepancy (MMD) based two-sample test (Gretton et al., ),
which compares two distributions based on their embeddings in some reproducing kernel Hilbert space (RKHS).
The estimated MMD serves as a measure of the extent to which the perturbations violate the disjointedness.
Interestingly, MMD test amounts to test Eq. 12 on a “most discriminative” feature map $h$ in the RKHS.
We provide a self-contained introduction about MMD and kernel mean embedding of
distributions in Appendix C.

Finally, it is helpful to consider a concrete generative process that yields disjoint perturbations.
Unlike the separability formulation, which assumes the existence of the Markov factorization of the latent distribution
(Assumption 3.4), the disjointedness models the latent distribution as a finite mixture.
In this framework, non-interacting perturbations intervene on different mixing components.

###### Assumption 3.8 .

The latent distribution $p_{Z}(z)$ admits the form of an $L$-component mixture:

$$ $\displaystyle p_{Z}(z|T_{1},\dots,T_{n})=\sum_{l=1}^{L}w_{l}\cdot p_{Z_{l}}(z| \operatorname{Pa}(z_{l})),\quad\text{where }(w_{1},\dots,w_{L})\in\Delta_{L}.$ (14) $$

In addition, perturbations do not intervene the mixing weights $(w_{1},\dots,w_{L})$.

###### Theorem 3.9 .

Suppose that Assumptions 3.1 and 3.8 hold. Then $\delta_{i},\delta_{j}$ are disjoint
if $\operatorname{Ch}(T_{i})\cap\operatorname{Ch}(T_{j})=\emptyset$.

###### Definition 3.7 .

###### Assumption 3.8 .

###### Theorem 3.9 .

### 4 Selecting perturbation pairs to efficiently discover interactions

With principled frameworks for testing separability and disjointedness in place, the natural next question is *how to select the experiments to run*? The space of possible pairs of perturbations is typically too large to perform experiments on all pairs as each experiment is costly. For example, pairwise knockouts of all 20 000 genes on the human genome would require approximately 200 million experiments (not including replicates). In this section we discuss how tools from experimental design and bandits can be used to efficiently select pairs of perturbations likely to have interactions.

###### Selection of perturbation pairs as active matrix completion.

The testing frameworks discussed in Section 3 prescribe test statistics which can be used to detect pairwise interactions. The test statistics, however, require samples from $p(x|\delta_{i,j})$, which entails running pairwise perturbation experiments. We are thus interested in developing an approach for selecting pairs of perturbations which are *likely* to have high values for the test statistics and are thus likely to reveal pairwise interactions.

The test statistics for pairs of perturbations can be viewed as an (unknown) symmetric matrix $\mathbf{R}\in\mathbb{R}^{n\times n}$, where each entry $\mathbf{R}_{i,j}$ contains the value of the test statistic that we will observe if we run perturbation $\delta_{ij}$. $\mathbf{R}$ is unknown a priori but we are allowed to select entries to observe. This can be viewed as an *active matrix completion* problem (Chakraborty et al., ). Active matrix completion is a variant of the standard matrix completion problem (Laurent, ) where values of entries of the matrix can be sequentially queried. Framing the problem of selecting perturbation pairs as an active matrix completion allows us to leverage existing efficient algorithms.

Figure: Algorithm 1 ASD for selecting perturbation pairs

###### Adaptive sampling for discovery.

In particular, we use the framework of *adaptive sampling for discovery* (ASD; Xu et al., ) which provides a bandit-based approach for active matrix completion. ASD involves using *information directed sampling* (IDS; Russo and Van Roy, ) in a “discovery” setting, where an action is only selected once. This is an instantiation of the general sleeping experts setting (Kanade et al., ), where the set of available actions shrinks every round.

Let $\Delta$ denote the action space of possible experiment designs, which in this case is the set of perturbations $\Delta:=\{\delta_{ij}:i,j\in n,i>j\}$, where $n$ is the total number of distinct perturbations.To simplify notation, we denote the perturbation pair $i,j$ selected at step $k$ as $a^{(k)}:=(i^{(k)},j^{(k)})\in\Delta$. $\mathcal{D}(\Delta)$ is the set of possible (categorical) distributions defined over $\Delta$. Each time an action, $a^{(k)}$, is selected, the corresponding element of the (unknown) reward matrix, $\mathbf{R}$, is revealed to the agent. In our setting this reward matrix corresponds to the test statistic from either Section 3.1 or 3.2 for each pair of perturbations. Let $H_{t}=((a^{(k)},\mathbf{R}_{i^{(k)},j^{(k)}}))_{k=1}^{t-1}$ denote the history of actions and their corresponding rewards until round $t$, and $\Delta_{t}$ denotes the set of remaining actions at round $t$; i.e. the pairs of perturbations that we have not yet tested experimentally. A policy $\pi$ is defined as a map from $H_{t}$ to $\mathcal{D}(\Delta_{t})$. The IDS policy, $\pi_{\text{IDS}}$ maintains a posterior distribution over $\mathbf{R}$ given the data observed up to round $t$, which we denote $p(\mathbf{R}\mid H_{t})$.

We can describe the sub-optimality of any action with respect to a set of beliefs by comparing the action’s reward to that of the best action that could have been selected at time $t$, under the agent’s current posterior over the reward matrix. Intuitively, we can evaluate this by sampling a plausible reward matrix from our posterior, $\dot{\mathbf{R}}\sim(\mathbf{R}|H_{t})$, and then comparing the reward from action, $a$, to the reward an agent would have received from selecting the optimal action, $a^{*}=\operatornamewithlimits{arg\,max}_{a\in\Delta_{t}}\dot{\mathbf{R}}(a)$; where $\dot{\mathbf{R}}(a^{(k)}):=\dot{\mathbf{R}}_{i^{(k)},j^{(k)}}$ for $a^{(k)}=(i^{(k)},j^{(k)})$. This is known as the expected instantaneous regret incurred by an action and is defined as,

$$ $\displaystyle\Delta_{t}(a)=\mathbb{E}_{\dot{\mathbf{R}}\sim p(\mathbf{R}|H_{t} )}[\dot{\mathbf{R}}(a^{*})-\dot{\mathbf{R}}(a)].$ (15) $$

Additionally, we define the information gain about the top $T-t+1$ remaining actions as follows:

$$ $\displaystyle g(a)=MI(a_{t,1}^{*},\dots,a_{t,T-t+1}^{*};\mathbf{R}(a)\mid H_{t },\delta_{t}=a).$ (16) $$

Algorithm 1 summarizes the algorithmic procedure. The algorithm operates over a series of $T$ rounds. In each round, the first step is to estimate the posterior distribution over $\mathbf{R}$ given the data observed thus far. Next step involves selecting a batch of perturbations based on the information ratio. The IDS policy at round $t$ can then be computed by minimizing the *information ratio* $\Psi$:

$$ $\displaystyle\pi_{\text{IDS}}\in\arg\min_{\pi\in{\cal D}(\Delta_{t})}\Psi_{\pi ,t}\coloneqq\frac{(\Delta_{t}^{\top}\pi)^{\lambda}}{g_{t}^{\top}\pi}$ (17) $$

where $\lambda$ is a parameter controls the tradeoff between lower instant regret (exploitation) and higher information gain (exploration). As $g_{t}(a)$ is intractable to compute in general, following Russo and Van Roy (); Xu et al. () we use an approximation, replacing $g_{t}$ with the conditional variance $v_{t}(a)=\text{Var}_{t}(\mathbb{E}[\mathbf{R}_{a}|a_{1}^{*},a])$. The conditional variance is a lower bound on the information gain $g_{t}(a)\geq v_{t}(a)$ and can thus replace the mutual information (since we are interested in maximizing it). After selecting the perturbation pair $i,j$, we compute the pairwise test statistic, $\mathbf{R}_{i,j}$ using data from the experimental outcomes. We add these test statistics to our reward matrix, update our posteriors and then continue on to the next round.

Without any structural assumptions on $\mathbf{R}$, linear regret is unavoidable for any bandit algorithm. We thus make an assumption that $\mathbf{R}$ is low-rank with a Gaussian prior on the columns. Under this assumption, ASD achieves sublinear regret (Xu et al., ).

###### Batching.

In high-throughput experimental screens, it is possible to run multiple experiments in parallel at the cost of a single experiment. So instead of a single action $\delta_{t}$, we can select a set of actions $\{\delta_{t}^{1},\dots,\delta_{t}^{b}\}$ where $b$ is the number of experiments we can run in parallel. However, as there are no efficient algorithms for combinatorial bandits in the discovery setting, we resort to a simple greedy scheme to select batches with ASD. Specifically, instead of picking a single action which minimizes the information ratio, we pick $b$ actions with the lowest information ratio.

### 5 Experiments

Our experiments aim to address three objectives: (1) verifying that our theoretical claims about
interactions are detectable on known synthetic tasks; (2) evaluating the test
statistics’ ability to recover known biological relationships from real
pairwise perturbation experiments; and (3) assessing our active learning
pipeline’s efficiency in detecting interactions.

In all the experiments,
our MMD-based tests used the RBF and Matern 2.5
kernels, and we chose bandwidth using median heuristics.
Unless otherwise stated, we estimated the KL score by first learning the three log-density ratios
$\log\frac{p(x|\delta_{ij})}{p(x|\delta_{0})},\log\frac{p(x|\delta_{i})}{p(x|
\delta_{0})},\log\frac{p(x|\delta_{j})}{p(x|\delta_{0})}$ using contrastive learning (Hermans et al., , NRE),
and then obtaining the KL estimates via
the smoothed mutual information “lower-bound” estimator (SMILE) (Song and Ermon, ) with clipping parameter $\tau=5$.
Detailed explanation about the NRE log-density ratio estimator and SMILE are provided in Appendix B.
Additional experimental details are provided in Appendix D.

#### 5.1 Testing on synthetic setting

We validated the separability test on two synthetic interventional examples:
one with 3-dimensional tabular data and one with images (of size $3\times 128\times 128$).
We also validated the disjointedness test on an interventional tabular example.
In each example, we generated observations by sampling from a latent distribution $p_{Z}$ that obeys a DAG structure,
followed by a mapping $g(\cdot)$ that maps the latent samples to the observations.
We then estimated the test statistics for the corresponding test for each pair of perturbations.
Detailed descriptions of the data generating processes for the three examples, the DAG structures of the latent distributions
$p_{Z}$ and the mapping from latents to observations $g(\cdot)$,
are provided in Section D.1.

Fig. 2 shows the synthetic results for the separability test on both
the tabular and image data. The results indicate that, in both examples, our estimated KL score
accurately characterized the separability relationships between perturbations:
inseparable pairs result in large KL scores, while separable
pairs result in small scores. The tests for disjointedness are shown in Fig. 3,
demonstrating that the MMD test accurately identified failures of
disjointedness and is relatively insensitive to the choice of kernel.

Figure: Figure 2: Separability testing on both the synthetic tabular data using KNN-based KL estimator (*left*) and NRE-based KL estimator (*middle*), and the synthetic images (*right*); brighter colors suggest stronger interactions. Ground truth interacting pairs for both examples are A-B and C-D, which are correctly identified.
Refer to caption: x2.png

Figure: Figure 3: Disjointedness testing on synthetic example using the MMD-based statistics with a Matern 2.5 kernel (*left*) and an RBF kernel (*right*); brighter colors suggest stronger interactions. The ground truth interacting pairs are D-E and F-G, which are correctly identified by the test.
Refer to caption: x5.png

#### 5.2 Testing on Biological interactions

Figure: Figure 4: *(left)* Pairwise separability scores between different CRISPR guides of two genes, TSC2 and MTOR. Missing pairs means that the data for corresponding pairwise combination was not collected in our biological experiment. Guides targeting the same gene show high KL scores (red), while guides targeting different genes show low scores (blue). *(right)* Random samples of the actual single cell images used in these experiments. Note that detecting the presence or absence of an interaction is extremely difficult, even for trained experts.
Refer to caption: x7.png

In order to evaluate whether our separability test can recover known biological interactions, we first ran the following test.
In gene knockout experiments, one can target the same gene with multiple different CRISPR guides, each of which cuts the gene in different places, but (at least in theory) results in the same gene being knocked out.
Intuitively, guides targeting the same gene should show high scores in the separability test because they are targeting the same latent variable, while guides targeting different genes should show lower scores (assuming the genes are on distinct pathways).
While two different guides that target the same gene are not usually run in a single well, we have in our dataset a couple of examples of this data for two genes (TSC2 and MTOR). For this experiment, we used single-cell cell painting images and tested the separability between pairs of guides.
The matrix of the separability scores displayed in Fig. 4 was consistent with what we would expect: we observed strong interaction scores between guides targeting the same gene (e.g. MTOR guide 3 with both MTOR guides 1 and 2), and much weaker scores for the interaction between the MTOR and TSC2 targeting guides. It is worth noting that both MTOR and TSC2 affect many systems within the cell, so we should expect some interaction, but the fact that the interaction was far smaller than the interaction score for guides targeting the same gene is very encouraging.

We then evaluated our testing approaches on a collection of 50 genes with the
goal of detecting gene interactions. The dataset was collected by performing
CRISPR knockouts on all pairs from a set of 50 selected genes, resulting in
1,225 gene pairs. The targeted genes have a bias towards known
gene-gene interactions.
We performed in vitro double gene knockout experiments on HUVEC cells using
three CRISPR guides per gene. We label perturbations according to the targeted gene, aggregating the effects of the individual CRISPR guides. The experimental protocol and data preprocessing followed the procedures
described by Fay et al. () and Sypetkowski et al. (),
respectively.

We computed the test statistics using 1024-dimensional embeddings of the
original cell painting images extracted from a pre-trained masked autoencoder
(He et al., ; Kraus et al., ).
Fig. 5 shows the matrices of the test statistics for all
perturbation pairs; each entry represents the pairwise interaction scores.
Qualitatively, both the disjointedness statistics and separability statistics effectively
uncover plausible biological relationships.
Many genes in the apoptosis pathway (programmed cell death, e.g., gene 2–3: BAX, BCL2L1) and the proteasome (protein degradation e.g., gene 28–30: PSMA1, PSMB2, PSMD1)
show high scores that are expected from synthetic lethal relationships. Synthetic lethality (SL) is a complex phenomenon in which simultaneous inactivation of specific gene combinations leads to cell death or extreme sickness, while individual perturbations have little effect. Apoptosis, controlling cell survival, is tightly regulated; disruption of the anti-apoptotic BCL-2 family (BCL2, BCL2L1, MCL1) (Kale et al., ) interferes with critical barriers against cell death. Proteasome function is similarly regulated, as cells must maintain a delicate balance of different proteins (Rousseau and Bertolotti, ). Apoptosis and proteasome members are often found in SL screens in cancer cells (Li et al., ; Ge et al., ; Han et al., ; Cron et al., ; Steckel et al., ; Das et al., ) because cancer genomes accumulate mutations that overcome weakened cellular buffering capabilities.

As displayed in Fig. 5, the MMD-based statistics show strong signals of interactions with proteasome components; the proteasome helps control global protein levels so we would expect to see it interacting with many different pathways, some essential.
The separability score provides less clear patterns but highlights several gene pairs that are known or expected to interact physically or genetically,
e.g., BCL2L1-MCL1 (gene 3-8) (Shang et al., ; Carter et al., ), BAX-BCL2L1 (gene 2-3) (Lindqvist and Vaux, ), BCL2L1-PSMD1 (gene 3-30) (Craxton et al., ), and PSMB2-PSMD1 (gene 29-30) (Voutsadakis, ).

Figure: Figure 5: Pairwise interaction scores from a selected set of $50$ genes using disjointedness test with Matern 2.5 kernel and RBF kernel (*left* and *middle*, respectively), and separability test (*right*); brighter colors suggest stronger interactions. Genes from the same pathway are ordered adjacently. The associated pathways of each selected gene are described in Table 1 of Section D.2.1.
Refer to caption: x8.png

#### 5.3 Interaction discovery with active matrix completion

Figure: Figure 6: Empirical results on the active pairwise experiments selection for the gene-gene interaction detection example using MMD based test statistics. The solid lines represent the mean performance, whereas the shaded region represent all the runs (min-max). IDS (ours) outperforms all the baselines significantly in terms of top scoring pairs discovered as well as the regret. In terms of known interactions, IDS still outperforms the baselines by a small margin.
Refer to caption: x11.png

Finally, we showcase how to use our test statistics to adaptively select the pairwise
experiment to run. As introduced in Section 4, the adaptive experimental design
is framed as an active matrix completion problem (Algorithm 1).
We applied Algorithm 1 on the gene-gene interaction detection example using the MMD-based scores.

###### Baselines.

To evaluate the effectiveness of our automated selection method, we compare against
selection with random policy, upper confidence bound
(UCB), Thompson sampling (TS) and uncertainty sampling (US). To instantiate UCB
in the discovery setting, we use the uncertainty from the low-rank matrix
posterior in place of the counts used in the standard multi-armed bandit
setting and mask the actions once they are selected. Similarly, in US we pick
pairs based solely on the uncertainty of the posterior. TS is instantiated in
the same way as IDS, but the pairs are selected to minimize only the instant
regret. Further details are discussed in Appendix D.

###### Evaluation metrics.

We evaluate each approach using three different metrics given a budget of 50
rounds, each with a batch size of 10 resulting in a total of 500 experiments
covering $~{}50\%$ of all possible pairs of the gene set.
First, we look at the fraction of the pairs with the top $5$ percentile of scores recovered by
the algorithm capturing the ability to explore the high scoring regions well.
Next, we also evaluate the regret of each algorithm with respect to an optimal
policy with access to the score matrix and acquires the highest scoring pairs
at each round. Finally, we evaluate the number of known biological relations
that appear in CORUM (Giurgiu et al., ), StringDB
(Szklarczyk et al., ), Signor (Licata et al., ) and hu.MAP
(Drew et al., ) to see how many each method is able to recover.

###### Results.

Fig. 6 (and Fig. 8 in Section D.3) illustrate the empirical results. We
observe that IDS discovers all the pairs with the top-5 percentile scores,
whereas all the baselines barely recover half of the top pairs. Even in terms
of the regret, IDS outperforms the baselines, with random performing the worst.
This demonstrates that IDS is able to exploit the low-rank structure in the
reward matrix effectively. Regarding known interactions discovered,
IDS and TS outperform other methods by a slight margin.
After 50 rounds, we observe that IDS and TS achieve between $12$ and $15\%$ improvement over the other baselines in the number of known biological interactions.
The significant difference between the performance on the fraction
of top pairs and the number of known relations prompts questions about the
correlation between the prediction errors and the biological interactions.

### 6 Discussion

This paper presented a method for efficiently detecting interactions between
perturbations. We present interactivity scores for both separability and disjointedness which both allow us to measure interactions between latent variables, and run
experiments only on pairs of perturbations for which we will fail to
compositionally generalize.
From an active learning perspective, disjointedness is powerful: when we can learn a good posterior over this test statistic, we effectively have a “confidence score” for whether two perturbations are likely to compose additively. This allows us to dramatically reduce the number of experiments we run by only experimenting where embeddings are unpredictable.
Separability gives a more intuitive notion of (in-)dependence in that we directly test for whether two perturbations interact in latent space. It will be interesting to explore to what extent this can be used to recover the target of interventions.

###### Limitations

While many known interactions scored highly according
to our test statistics, the overall correlation between known
interactions and this metric was relatively low. It is not clear whether this
is the result of a lack of specificity, or whether this is because we are in
fact discovering real relationships that are not known to biology. To test this
would require additional experimentation with orthogonal assays.
Additionally, the results for the separability tests depended on the quality of the KL estimator. We used the SMILE estimator which is somewhat sensitive to the choice of the clipping parameter, $\tau$. In our experiments on images, a default of $\tau=5$ seemed to work well, but if we want this test to be applicable across a variety of modalities, we need a more robust method for choosing this hyper-parameter.

### 7 Acknowledgements

We are extremely grateful for the many discussions with colleagues at Recursion and Valence Labs, and external collaborators that lead to this work.
The original design of the benchmark dataset was designed by Marta Fay and the experiments were managed and run by Jordan Finnell, Brandon Mendivil, Kate Brown and Vicky Chen. Thank you to Nathan Lazar for his help in preparing the data and helpful discussions throughout the project.
The pairwise interaction testing section was strongly influenced by discussions with everyone at the Bellairs Causal Workshop; special thanks to Victor Veitch, Kartik Ahuja and Yixin Wang for their input.
Finally, thank you to Cian Eastwood, Johnny Xi and Jana Osea for the helpful discussions and feedback on the work.

## Appendix

### Appendix A Proofs

###### Proof of Theorem 3.6 .

By Assumptions 3.4 and 3.2, we can apply change of variable formula to obtain the
log density of $p_{X}(x|T_{1},\dots,T_{n})$:

$$ $\displaystyle\log p(x|T_{1},\dots,T_{n})=\sum_{l=1}^{L}\log p_{z_{l}}([g^{-1}( x)]_{l}|\operatorname{Pa}(z_{l}))+\log|\det\nabla g^{-1}(x))|.$ (18) $$

Here $[\cdot]_{l}$ denote the projection operator that maps $z\in{\mathcal{Z}}$ to the subspace on which $z_{l}$ lies, i.e.,
$\forall z\in{\mathcal{Z}},[z]_{l}=z_{l}$.
Recall that we interpret $z$ as the concatenation of all latent factors $(z_{1},\dots,z_{L})$.
If $T_{i}$, $T_{j}$ are causally independent, perturbation $\delta_{i},\delta_{j}$
will intervene different terms in the above summand. Without lose of generality,
suppose $\delta_{i}$ and $\delta_{j}$ intervene $z_{l_{i}}$ and $z_{l_{j}}$ repsectively.
Then,

$$ $\displaystyle\bar{\ell}_{X}(x|\delta_{i})$ $\displaystyle=\log p(x|\delta_{i})-\log p(x|\delta_{0})$ (19) $\displaystyle=\log p_{z_{l_{i}}}([g^{-1}(x)]_{l_{i}}|\delta_{i})-\log p_{z_{l_ {i}}}([g^{-1}(x)]_{l_{i}}|\delta_{0}).$ (20) $$

where the equality follows from the modularity assumption that the intervention only affects $l_{i}$.

Similarly, we obtain that

$$ $\displaystyle\bar{\ell}_{X}(x|\delta_{j})$ $\displaystyle=\log p_{z_{l_{j}}}([g^{-1}(x)]_{l_{j}}|\delta_{j})-\log p_{z_{l_ {j}}}([g^{-1}(x)]_{l_{j}}|\delta_{0}),$ (21) $\displaystyle\bar{\ell}_{X}(x|\delta_{ij})$ $\displaystyle=\log p_{z_{l_{i}}}([g^{-1}(x)]_{l_{i}}|\delta_{i})-\log p_{z_{l_ {i}}}([g^{-1}(x)]_{l_{i}}|\delta_{0})$ (22) $\displaystyle+\log p_{z_{l_{j}}}([g^{-1}(x)]_{l_{j}}|\delta_{j})-\log p_{z_{l_ {j}}}([g^{-1}(x)]_{l_{j}}|\delta_{0}),$ (23) $$

which completes the proof.
∎

###### Proof of Theorem 3.9 .

Provided with Assumption 3.1, $\delta_{i},\delta_{j}$ being disjoint (Definition 3.7) is equivalent to
the same additivity in the intervened latent distributions, i.e.,

$$ $\displaystyle\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{ eq:latentexclusive}}{e}q:latentexclusive}p_{Z}(z|\delta_{ij})-p_{Z}(z|\delta_{ 0})=p_{Z}(z|\delta_{i})-p_{Z}(z|\delta_{0})+p_{Z}(z|\delta_{j})-p_{Z}(z|\delta _{0}).$ (24) $$

Under the mixture model assumption Assumption 3.8, if $T_{i}$, $T_{j}$ are causally independent,
if $\operatorname{Ch}(T_{i})\cap\operatorname{Ch}(T_{j})=\emptyset$,
perturbation $\delta_{i},\delta_{j}$ will intervene distinct mixing components.
Without lose of generality, suppose $\delta_{i}$ and $\delta_{j}$ intervene $p_{Z_{l_{i}}}$ and $p_{Z_{l_{j}}}$ repsectively.

$$ $\displaystyle p_{Z}(z|\delta_{i})-p_{Z}(z|\delta_{0})$ $\displaystyle=w_{l_{i}}\cdot(p_{Z_{l_{i}}}(z|\delta_{i})-p_{Z_{l_{i}}}(z| \delta_{0})),$ (25) $\displaystyle p_{Z}(z|\delta_{j})-p_{Z}(z|\delta_{0})$ $\displaystyle=w_{l_{j}}\cdot(p_{Z_{l_{j}}}(z|\delta_{j})-p_{Z_{l_{j}}}(z| \delta_{0})),$ (26) $\displaystyle p_{Z}(z|\delta_{ij})-p_{Z}(z|\delta_{0})$ $\displaystyle=w_{l_{i}}\cdot(p_{Z_{l_{i}}}(z|\delta_{i})-p_{Z_{l_{i}}}(z| \delta_{0}))+w_{l_{j}}\cdot(p_{Z_{l_{j}}}(z|\delta_{j})-p_{Z_{l_{j}}}(z|\delta _{0})),$ (27) $$

which shows Eq. 24 and hence completes the proof.
∎

###### Proof of Theorem 3.6 .

###### Proof of Theorem 3.9 .

### Appendix B Sample-based estimation of KL divergence

Sample-based estimator of KL-divergence is a challenging task, particularly for high dimensional observations.
In this section, we present the estimation procedure we used in our experiments.
As mentioned in Section 5, to estimate the KL-divergence,

$$ $\displaystyle\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{ eq:kl}}{e}q:kl}\mathrm{D_{KL}}(p||p_{i})=\mathbb{E}\left[\log\frac{p(X|\delta_ {0})}{p(X|\delta_{i})}\right],\quad X\sim p(X|\delta_{0}),$ (28) $$

we used a two-step procedure:

- 1.
We first estimate the log-density ratio
$\displaystyle\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{
eq:ldr}}{e}q:ldr}\log\frac{p(x|\delta_{0})}{p(x|\delta_{i})}$
(29)
via a neural ratio estimator (NRE) based on a contrastive learning objective [Hermans et al., , NRE].
- 2.
Then, instead of taking naive Monte Carlo estimates of the KL-divergence, we adopt the SMILE estimator [Song and Ermon, ] using the learned log-density ratio.

Sections B.1 and B.2 describe the two steps respectively.

#### B.1 Contrastive neural ratio estimator

There are various methods to learn the log-density ratio between a single pair
of distributions. For example, the optimal discriminator between two
distributions is related to their density ratios [Goodfellow et al., , Proposition
1]. However, a common structure in our applications
involves a large set of perturbations, and training a classifier for each
perturbation class (against the control group) would be very cumbersome and not
data efficient. Therefore, we consider a contrastive learning model
[Hermans et al., , NRE], which trains a binary classifier to
distinguish the joint data distribution
$p(x,c)$ from the product of the marginals
$p(x)p(c)$, where $c$ denotes the perturbation class.

The training objective is as follows:

$$ $\displaystyle\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{ eq:nre}}{e}q:nre}\theta,w\in\arg\min_{\theta,w}-\frac{1}{2B}\left[\sum_{b=1}^{ B}\log(1-\text{Sigmoid}(f_{\theta,w}(x^{(b)},c^{(b)}))+\sum_{b^{\prime}=1}^{B} \log(\text{Sigmoid}(f_{\theta,w}(x^{(b^{\prime})},c^{(b^{\prime})}))\right],$ (30) $$

where $x^{(b)},c^{(b)}\sim p(x)p(c)$ and $x^{(b^{\prime})},c^{(b^{\prime})}\sim p(x,c)$, and

$$ $\displaystyle f_{\theta,W}(x,c)=\text{Encoder}_{\theta}(x)^{T}W_{c}.$ (31) $$

Given infinite training samples and flexibility of the neural network $f$, the optimal $f^{\star}(x,c)=\log\frac{p(x,c)}{p(x)p(c)}=\log\frac{p(x|c)}{p(x)}$.
And we can then obtain Eq. 29 via $f^{\star}(x,\delta_{i})-f^{\star}(x,\delta_{0})$.

In our examples, we estimate our log-density ratios by training the NRE objective on all perturbation classes.

#### B.2 Smoothed Mutual Information “Lower-bound” Estimator (SMILE)

After obtaining the log-density ratio estimator for Eq. 29, there are several options for estimating the KL-divergence. We provide a short review here on the various strategies [Ghimire et al., , Belghazi et al., , Hjelm et al., , Song and Ermon, ], and explain why we opt for the SMILE estimator [Song and Ermon, ].

The most straightforward estimates of the KL-divergence Eq. 28 is by the Monte
Carlo estimates based on samples from $p(X|\delta_{0})$, i.e.,

$$ $\displaystyle\mathrm{D_{KL}}(P(X|\delta_{i})||P(X|\delta_{0}))\approx\frac{1}{ N}\sum_{i=1}^{N}f(X_{i}),\quad\text{where }f(\cdot)\text{is the estimated}\log \frac{p(\cdot|\delta_{i})}{p(\cdot|\delta_{0})}.$ (32) $$

However, it is noted that the variance of this estimator is often huge
[Song and Ermon, , Ghimire et al., ], making the estimated KL
unreliable in practice.

Belghazi et al. [] proposed to estimate the KL-divergence based on its Donsker-Varadhan representation [Donsker and Varadhan, ], i.e.,

$$ $\displaystyle\mathrm{D_{KL}}(p||q)=\sup_{f:\Omega\to\mathbb{R}}\mathbb{E}_{p}[ f]-\log\mathbb{E}_{q}[\exp(f)],$ (33) $$

where the supremum is taken over all functions $f$ such that the two expectations are finite.
Notice that the optimal $f$ is indeed achieved as the log-density ratio between $p$ and $q$.
In practice, one can parameterize $f$ using some neural network and maximize the above objective to approximate the KL [Belghazi et al., ]. However, the stochastic gradient estimator of the above objective is generally biased [Belghazi et al., ], making the optimization less stable.
Therefore, it is also recommended to
first learn the log-density ratio $\log\frac{p}{q}$ as $f$, and then estimate the KL-divergence [Hjelm et al., , Song and Ermon, ] via

$$ $\displaystyle\text{KL}_{\text{MILE}}(f)=\mathbb{E}_{p}[f]-\log\mathbb{E}_{q}[ \exp(f)],$ (34) $$

which technically gives a lower-bound of the KL divergence if the log-density ratio is not well-estimated. We would refer this estimator as the mutual information lower-bound estimator (MILE).
MILE is often shown to have better performance than the naive Monte Carlo estimator [Belghazi et al., , Hjelm et al., , Song and Ermon, ].

However, MILE also suffers from high variance issues [Song and Ermon, ] particularly when the learned $f$ has large values in the tail of $q$.
Song and Ermon [] proposed a smoothed version of MILE, named as SMILE, by clipping the learned log-density ratios $f$ between $-\tau$ and $\tau$, i.e.,

$$ $\displaystyle\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{ eq:smile}}{e}q:smile}\text{KL}_{\text{SMILE}}(f,\tau)=\mathbb{E}_{p}[f]-\log \mathbb{E}_{q}[\text{clip}(\exp(f),\exp(-\tau),\exp(\tau))].$ (35) $$

$\text{KL}_{\text{SMILE}}(f,\tau)$ converges to $\text{KL}_{\text{MILE}}(f)$ as $\tau\to\infty$, but smaller $\tau$ significantly reduces the variance of the MILE estimator.

In our experiments, we use the SMILE estimator with the clipping parameter $\tau$ set to be $5$.

### Appendix C MMD and kernel mean embedding

#### C.1 Kernel mean embedding

We provide here a minimal overview of the kernel mean embedding of probability measures [Muandet et al., ].
Given a feature map $\phi:{\mathcal{X}}\to{\mathcal{F}}_{\phi}$, where ${\mathcal{F}}_{\phi}$ is some Hilbert space (sometimes called the feature space).
This feature map $\phi$ defines a kernel:

$$ $\displaystyle k_{\phi}:{\mathcal{X}}\times{\mathcal{X}}\to\mathbb{R},\qquad k_ {\phi}(x,y)=\langle\phi(x),\phi(y)\rangle,$ (36) $$

where $\langle\cdot,\cdot\rangle$ denotes the inner product of ${\mathcal{F}}_{\phi}$.
Kernel $k_{\phi}$ defined this way induces a space of functions—${\mathcal{H}}_{\phi}$—from
${\mathcal{X}}$ to $\mathbb{R}$, which is a *reproducing kernel Hilbert space* (RKHS).
The name RKHS comes from a special property of ${\mathcal{H}}_{\phi}$, called the *producing property*:

$$ $\displaystyle\forall f\in{\mathcal{H}}_{\phi},\qquad\langle f,\phi(x)\rangle_{ {\mathcal{H}}_{\phi}}=f(x).$ (37) $$

We abuse the notation here by interpreting $\phi(x)$ as a function in
${\mathcal{H}}_{\phi}$ instead of a real value.
A special instance of the reproducing property is that

$$ $\displaystyle\langle\phi(x),\phi(y)\rangle=k_{\phi}(x,y).$ (38) $$

Precisely, we are using the following representation of $\phi(x)$:

$$ $\displaystyle\phi(x):{\mathcal{X}}\to{\mathcal{H}}_{\phi},\qquad\phi(x)(\cdot) =k_{\phi}(x,\cdot)\in{\mathcal{H}}_{\phi}.$ (39) $$

Here $k_{\phi}(x,\cdot)\in{\mathcal{H}}_{\phi}$ is guaranteed by the definition of ${\mathcal{H}}_{\phi}$.

In what follows, we use ${\mathcal{M}}^{1}_{+}({\mathcal{X}})$ to denote the space of probability measures over ${\mathcal{X}}$.
We can embed any probability measure to the RKHS.

###### Definition C.1 (Kernel mean embedding of probability measure) .

The *kernel mean embedding* of a probability measure $\mathbb{P}\in{\mathcal{H}}_{\phi}$ is defined via the following mapping:

$$ $\displaystyle\mu_{\phi}:{\mathcal{M}}^{1}_{+}({\mathcal{X}})\to{\mathcal{H}}_{ \phi},\qquad\mu_{\phi}(\mathbb{P})(\cdot)=\int_{{\mathcal{X}}}\phi(x)(\cdot) \mathbb{P}(\mathrm{d}x)=\mathbb{E}_{X\sim\mathbb{P}}\left[\phi(X)\right](\cdot).$ (40) $$

We denote the kernel mean embedding of $\mathbb{P}$ by $\mu_{\phi}(\mathbb{P})$.

###### Proposition C.2 (c.f Eq. (3.29) of Muandet et al. [ 2017 ] ) .

For all $\mathbb{P},\mathbb{Q}\in{\mathcal{M}}^{1}_{+}({\mathcal{X}})$,

$$ $\displaystyle\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{ eq:mmdl2norm}}{e}q:mmdl2norm}\operatorname{MMD}_{k_{\phi}}(\mathbb{P},\mathbb{ Q})=\left\|\mathbb{E}_{X\sim\mathbb{P}}\left[\phi(X)\right]-\mathbb{E}_{Y\sim \mathbb{Q}}\left[\phi(Y)\right]\right\|=\|\mu_{\phi}(\mathbb{P})-\mu_{\phi}( \mathbb{Q})\|.$ (41) $$

###### Proof of Proposition C.2 .

The second equality of Eq. 41 follows directly from the mapping defined in Definition C.1.
We then focus on proving the first equality.

$$ $\displaystyle\operatorname{MMD}_{k_{\phi}}(\mathbb{M}_{1},\mathbb{M}_{2})$ (42) $\displaystyle=\sup_{f\in{\mathcal{H}}_{\phi}:\|f\|\leq 1}\mathbb{E}_{X\sim \mathbb{P}}[f(X)]-\mathbb{E}_{Y\sim\mathbb{Q}}[f(Y)]$ (43) $\displaystyle=\sup_{f\in{\mathcal{H}}_{\phi}:\lVert f\rVert\leq 1}\mathbb{E}_{ X\sim\mathbb{P}}[\langle f,\phi(X)\rangle]-\mathbb{E}_{Y\sim\mathbb{Q}}[ \langle f,\phi(Y)\rangle]\quad\text{(reproducing property)}$ (44) $\displaystyle=\sup_{f\in{\mathcal{H}}_{\phi}:\lVert f\rVert\leq 1}\langle f, \mathbb{E}_{X\sim\mathbb{P}}[\phi(X)]\rangle-\langle f,\mathbb{E}_{Y\sim \mathbb{Q}}[\phi(Y)]\rangle\quad\text{(linearity of inner product)}$ (45) $\displaystyle=\sup_{f\in{\mathcal{H}}_{\phi}:\lVert f\rVert\leq 1}\langle f, \mathbb{E}_{X\sim\mathbb{P}}[\phi(X)]-\mathbb{E}_{Y\sim\mathbb{Q}}[\phi(Y)] \rangle\quad\text{(linearity of inner product)}$ (46) $$

Then the proof is completed by the observation that the inner product is maximized when

$$ $\displaystyle f=\left\{\mathbb{E}_{X\sim\mathbb{P}}[\phi(X)]-\mathbb{E}_{Y\sim \mathbb{Q}}[\phi(Y)]\right\}/\left\|\mathbb{E}_{X\sim\mathbb{P}}[\phi(X)]- \mathbb{E}_{Y\sim\mathbb{Q}}[\phi(Y)]\right\|.$ (48) $$

∎

Proposition C.2 essentially says that MMD between two distributions
is indeed the distance of mean embeddings of features.
It also tells us that $\operatorname{MMD}_{k_{\phi}}(\mathbb{P},\mathbb{Q})=0$ if and only if $\mu_{\phi}(\mathbb{P})=\mu_{\phi}(\mathbb{Q})$.
To be able to separate any two distributions via MMD, the kernel mean embedding must be an injective map,
in which case the feature map induces a *characteristic kernel*.

###### Definition C.3 (Characteristic kernel) .

$k$ is said to be characteristic on ${\mathcal{M}}^{1}_{+}({\mathcal{X}})$ if the kernel mean embedding,

$$ $\displaystyle\mu_{k}:{\mathcal{M}}^{1}_{+}({\mathcal{X}})\to{\mathcal{H}}_{k}, \qquad\mu_{k}(\mathbb{P})(\cdot)=\int_{{\mathcal{X}}}\phi(x)(\cdot)\mathbb{P}( \mathrm{d}x),$ (49) $$

is injective. In other words, $k$ is a characteristic kernel if and only if

$$ $\displaystyle\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{ eq:characteristickernel}}{e}q:characteristickernel}\mu_{k}(\mathbb{P},\mathbb{ Q})=0\Longleftrightarrow\mu_{k}(\mathbb{P})=\mu_{k}(\mathbb{Q}) \Longleftrightarrow\mathbb{P}=\mathbb{Q}.$ (50) $$

Here we change the subscript of $\mu$, ${\mathcal{H}}$ from the feature map to the kernel,
because from now on we do not necessarily work on explicit choice of feature maps;
many characteristic kernels do not have tractable feature maps (and mostly infinite dimensional).
Eq. 50 states that MMD is a metric on ${\mathcal{M}}^{1}_{+}({\mathcal{X}})$ if a characteristic kernel is used; otherwise, distinct distributions with the same kernel mean embedding cannot be separated by MMD.
Examples of characteristic kernels on ${\mathcal{M}}^{1}_{+}(\mathbb{R}^{d})$ are Gaussian kernels, Laplacian kernels, and the family of Matérn kernels.
We refer readers to Sriperumbudur et al. [] for a comprehensive
survey of the characteristic kernels.

Provided with a specified kernel, and samples from $\mathbb{P},\mathbb{Q}$, one can obtain an unbiased estimate of the squared population MMD;
see Gretton et al. [, Eq. (3)] for the detailed expression of the estimator.
An interesting property of the estimates of MMD is that the convergence is
dimension independent [Gretton et al., , Theorem 7],
although the penalization of the dimension in the context of kernel two sample
test lies in the reduction of power [Ramdas et al., ].

###### Definition C.1 (Kernel mean embedding of probability measure) .

###### Proposition C.2 (c.f Eq. (3.29) of Muandet et al. [ 2017 ] ) .

###### Proof of Proposition C.2 .

###### Definition C.3 (Characteristic kernel) .

#### C.2 Compare distributions via a fixed feature map h ℎ h italic_h

In some cases, one can compare $\mathbb{P},\mathbb{Q}$ by assessing on a fixed test function $h$.
For example, the L2 norm between the expectation of $h$, i.e.,

$$ $\displaystyle\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{ eq:l2rewardh}}{e}q:l2rewardh}\|\mathbb{E}_{X\sim\mathbb{P}}[h(X)]-\mathbb{E}_{ Y\sim\mathbb{Q}}[h(Y)]\|_{2}$ (51) $$

can be a crude measure on how different $\mathbb{P},\mathbb{Q}$. According to Proposition C.2,

$$ $\displaystyle\|\mathbb{E}_{X\sim\mathbb{P}}[h(X)]-\mathbb{E}_{Y\sim\mathbb{Q}} [h(Y)]\|_{2}=\mathrm{MMD}_{k_{h}}(\mathbb{P},\mathbb{Q}).$ (52) $$

The obvious limitation with a fixed choice of $h$ is that it’s typically not characteristic for all distributions, i.e.,

$$ $\displaystyle\|\mathbb{E}_{X\sim\mathbb{P}}[h(X)]-\mathbb{E}_{Y\sim\mathbb{Q}} [h(Y)]\|_{2}=0\mathchoice{\mathrel{\hbox to0.0pt{\kern 3.75pt\kern-5.27776pt$ \displaystyle\not$\hss}{\implies}}}{\mathrel{\hbox to0.0pt{\kern 3.75pt\kern-5 .27776pt$\textstyle\not$\hss}{\implies}}}{\mathrel{\hbox to0.0pt{\kern 2.625pt \kern-4.11108pt$\scriptstyle\not$\hss}{\implies}}}{\mathrel{\hbox to0.0pt{ \kern 1.875pt\kern-3.3333pt$\scriptscriptstyle\not$\hss}{\implies}}}\mathbb{P} =\mathbb{Q}.$ (53) $$

However, in practice, if the choice of $h$ is sufficient to identify discriminate the
collection of distributions that we care about, it is convenient to just examine Eq. 51.

In our specific application, interaction detection, we find the embedding vector that are
constructed from the final hidden layer of a classifier works reliably.
We assume that this classifier is trained optimally such that, $\mathbb{P}(\delta_{i}|X)=\sigma(w_{i}^{\top}h(X))$, where $\sigma:=\frac{\exp(x_{i})}{\sum_{j}\exp(x_{j})}$ is the softmax function, and that there are
sufficiently diverse labels to ensure that this representation is identified up
to a linear transformation; see Roeder et al. [] for details. Given a trained
classifier,
$\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar
\macc@nested@a 111{}_{i}:=\mathbb{E}[h(X)|\delta_{i}]-\mathbb{E}[h(X)|\delta_{
0}]$
is then the average of the embeddings associated with a particular knockout centered around the control wells.
We find in our empirical experiments that the metric $\|\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar
\macc@nested@a 111{}_{ij}-\macc@depth\char 1\relax\frozen@everymath{
\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{i}-\macc@depth\char 1
\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}_{j}
\|_{2}$
works well to identify highly interacting gene pairs.

We have not obtained theoretical arguments on justifying the use of optimal discriminator as the choice of test functions;
we left this for future development.

### Appendix D Experiment details

For all separability tests, we trained the NRE model on all perturbation
classes and then obtained the log-density ratios for each pair as introduced in Section B.1.
For the NRE training, we evaluated multiple model architectures and optimizer
step sizes, selecting the hyperparameter combination that yielded the best
training accuracy. The selected model architectures and optimizer step sizes
for each example are reported in the corresponding sections. We checkpointed
the best model at the optimal training accuracy for further inference to obtain
density ratio estimates and KL estimates. To optimize our models, we used ADAM [Kingma and Ba, ] with default
hyperparameter settings.

All our experiments were run on NVIDIA H100 GPUs.

#### D.1 Synthetic examples

For all the synthetic examples, we generated 20,000 i.i.d. samples for each
perturbation class, including both single and double perturbations. The
detailed data generation process for each example is provided below.

##### D.1.1 Synthetic tabular for separability test

The separability scores in this example are computed using two different
KL-divergence estimators: the simple KNN-based estimator and the NRE-based
estimation procedure described in Appendix B. We used a 3-layer MLP with
ReLU activation (hidden dimensions of 128 and 64) as the encoder for the NRE
density ratio estimator. The NRE model was trained using the ADAM optimizer
with a step size of 0.005 for 500 epochs and a batch size of 1024.

###### Latent distribution

The latent variable $Z$ consists of 3 independent one-dimensional variables $P_{1},P_{2},P_{3}$, i.e.,

$$ $\displaystyle P_{Z}(z_{1},z_{2},z_{3})=P_{1}(z_{1})\cdot P_{2}(z_{2})\cdot P_{ 3}(z_{3}),$ (54) $$

for which we consider 4 single perturbations in total, labelled as $A,B,C,D$ respectively, resulting $6$ pairwise perturbations.
Perturbations are applied by changing the distribution of one or more of latent variables.
In our setting, perturbations are classified as separable or inseparable.
Double perturbations of separable ones will intervene in two distinct latent
variables, while double perturbations of inseparable ones will intervene in the same latent variable.
For each node, the control (unperturbed) distribution is $\mathcal{N}(0,1)$, while the perturbed one—whether it’s single perturbed or doubly perturbed—becomes $\mathcal{N}(3,1)$.

The association between the perturbation class and corresponding intervened latent variable(s) is described as follows:

$$ $\displaystyle\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{ eq:syncdag}}{e}q:syncdag}\begin{aligned} P_{1}&\sim\begin{cases}\text{Normal}( 0,1)&\text{if unperturbed}\\ \text{Normal}(3,1)&\text{if perturbed by (A) or/and (B)}\end{cases}\\ P_{2}&\sim\begin{cases}\text{Normal}(0,1)&\text{if unperturbed}\\ \text{Normal}(3,1)&\text{if perturbed by (B)}\end{cases}\\ P_{3}&\sim\begin{cases}\text{Normal}(0,1)&\text{if unperturbed}\\ \text{Normal}(3,1)&\text{if perturbed by (C) or/and (D)}\end{cases}\end{aligned}$ (55) $$

As per Eq. 55, the only two inseparable pairs are A-B (both intervening $P_{1}$) and C-D (both intervening $P_{3}$).

###### Sampling process

To generate observations for each perturbation class, we first generate latent
samples from $P_{Z}$ based on the Eq. 55, and then transform the
latent samples via a diffeomorphism $g(\cdot)$. In this example, we chose
$g(\cdot)$ as a randomly initialized 7-layer Multi layer perceptron (MLP) with
LeakyReLU activations.

##### D.1.2 Synthetic image example for separability test

The synthetic images consist of three objects (three small colored balls) in
different locations and with various backgrounds. The positions of the objects
are encoded in a 3-dimensional latent variable $Z$, which follows the
identical DAG structure described in Eq. 55. Each coordinate of
$Z$ corresponds to the location distribution of an object, determining the
distribution of both the $x$ and $y$ coordinates of the object. Thus, a
perturbation affecting the location distribution of one object will intervene
in the distribution of both coordinates of that object. Background distortion
is controlled by random noise. Scenes are generated using a rendering engine
from PyGame, denoted as $g(\cdot)$. Example images are provided in
Fig. 7.

For the NRE model, we used a 5-layer convolutional network (with channels 32,
54, 128, 256, and 512) featuring batch normalization and leaky ReLU activations
as the encoder. The model was trained using the ADAM optimizer with a step size
of 0.0002 for 200 epochs and a batch size of 1024.

Figure: Figure 7: Example images of the interventional image data. 3 small balls with blue, purple, and red colors respectively are the targeted objects, whose locations are intervened by perturbations. Large balls with green and yellow colors form the background.
Refer to caption: extracted/5847700/figs/sync/example1.png

##### D.1.3 Synthetic tabular example for disjointedness test

The latent distribution is set to be a 6-component mixture distribution,

$$ $\displaystyle P_{Z}(z)=\frac{1}{8}P_{0}(z)+\frac{1}{8}P_{1}(z)+\frac{1}{8}P_{2 }(z)+\frac{1}{8}P_{3}(z)+\frac{1}{4}P_{4}(z)+\frac{1}{4}P_{5}(z).$ (56) $$

We define 7 single perturbations, labeled as A, B, C, D, E, F, G, respectively.
The control distributions (unperturbed distribution) of all mixture component are $\mathcal{N}(0,1)$.
Perturbations are applied by changing the distribution of one or more of the mixture components.
E.g., if perturbation $D$ is applied, $P_{4}\sim\text{Normal}(0,5)$.
In our setting, perturbations are classified as independent ones and interacting ones,
where double perturbation of independent ones will intervene two distinct mixture components,
while double perturbation of interacting ones will intervene the same mixture component.
The association of the perturbations and the corresponding component gets intervened is provided as follows:

The independent group consists of the mixture components $P_{0},P_{1},P_{2},P_{3}$:

$$ $\displaystyle\begin{aligned} P_{0}&\sim\text{Normal}(0,1)&\text{(not intervened)}\\ P_{1}&\sim\text{Normal}(5,1)&\text{(A)}\\ P_{2}&\sim\text{Normal}(10,1)&\text{(B)}\\ P_{3}&\sim\text{Normal}(-5,5)&\text{(C)}\end{aligned}$ (57) $$

The interacting group consists of the variables $P_{4},P_{5}$:

$$ $\displaystyle\begin{aligned} P_{4}&\sim\begin{cases}\text{Normal}(0,5)&\text{( D)}\\ \text{Normal}(-10,5)&\text{(E)}\\ \text{Normal}(10,5)&\text{(DE)}\end{cases}\\ P_{5}&\sim\begin{cases}\text{Cauchy}(15,1)&\text{(F)}\\ \text{Cauchy}(-15,1)&\text{(G)}\\ \text{Cauchy}(20,1)&\text{(FG)}\end{cases}\end{aligned}$ (58) $$

###### Sampling process

To generate latent samples, we draw from the specified distributions and the constructed mixture model. Let $n$ be the number of samples and $d$ the dimension. The samples are generated as follows:

$$ $\displaystyle\begin{aligned} &\text{if no intervention:}\quad\mathbf{X}\sim P_ {0}\\ &\text{if intervention is specified:}\quad(\mathbf{P}_{s},\mathbf{w}_{s})= \text{build_mixture}(\text{intervene})\\ &\text{draw categorical samples:}\quad\mathbf{z}\sim\text{Categorical}(\mathbf {w}_{s})\\ &\text{generate samples from each component:}\quad\mathbf{X}_{i}\sim\mathbf{P} _{s}[i]\end{aligned}$ (59) $$

We generated latents for all 7 single perturbations, and all double perturbations, and
then obtained the observations, which are used to perform tests on, by mapping the latent samples through a deterministic function $g(\cdot)$.
We chose $g(\cdot)$ as a randomly initialized 10-layer MLP with LeakyReLU activations.

#### D.2 Real data examples

##### D.2.1 Gene-gene interactions

For the separability test, we used a 3-layer MLP with ReLU activation (hidden
dimensions of 2048 and 256) as the encoder for the NRE density ratio estimator.
We trained the NRE model using the ADAM optimizer with a step size of $0.0001$
for 2500 epochs and a batch size of 16384.

Table 1 describes the corresponding pathways for each selected gene.

**Table 1: List of gene indices and their associated pathways**
| Gene Index | Pathway |
| --- | --- |
| gene0 | Amino acid sensing (mTOR pathway) |
| gene1 | Apoptosis |
| gene2 |  |
| gene3 |  |
| gene4 |  |
| gene5 |  |
| gene6 |  |
| gene7 |  |
| gene8 |  |
| gene9 |  |
| gene10 | Autophagy |
| gene11 |  |
| gene12 |  |
| gene13 |  |
| gene14 |  |
| gene15 |  |
| gene16 |  |
| gene17 |  |
| gene18 |  |
| gene19 |  |
| gene20 | ERAD (protein folding) |
| gene21 |  |
| gene22 |  |
| gene23 | Integrated Stress Response |
| gene24 | Microtubule |
| gene25 |  |
| gene26 |  |
| gene27 | PI3K-Akt signaling |
| gene28 | Proteasome |
| gene29 |  |
| gene30 |  |
| gene31 | Protein translation |
| gene32 | Protein translation (mTOR pathway) |
| gene33 | Ribosome |
| gene34 |  |
| gene35 |  |
| gene36 | Transcriptional regulation |
| gene37 | UPR (protein folding) |
| gene38 |  |
| gene39 |  |
| gene40 |  |
| gene41 |  |
| gene42 |  |
| gene43 | mTOR signaling |
| gene44 |  |
| gene45 |  |
| gene46 |  |
| gene47 |  |
| gene48 |  |
| gene49 | p53 signaling |

##### D.2.2 Guide-guide interactions

The single-cell painting images are derived from multi-cell images, with each
single-cell nucleus centered within a $32\times 32$ pixel box. The encoder
of the NRE model maps single-cell images of shape $(6,32,32)$ into a
128-dimensional feature vector. It consists of three convolutional blocks, each
comprising a Conv2D layer with a 3x3 kernel, BatchNorm2D, ReLU activation, and
MaxPool2D, progressively increasing the number of channels from 6 to 32, 64,
and 128 while halving the spatial dimensions at each max-pooling step. After
the convolutional layers, the output tensor of shape $(128,4,4)$ is
flattened to $(2048)$ and passed through two fully connected blocks, each
with a Linear layer, ReLU activation, and Dropout (with a dropout rate of 0.3),
transforming the feature size from 2048 to 256 and finally to 128. We train the
NRE model with the ADAM optimizer, using a step size of $0.00005$ for 5000
epochs and a batch size of 2048.

#### D.3 Active learning

To obtain samples from the posterior distribution over the low-rank reward
matrix (with rank $m$), we use stochastic variational
inference [Wingate and Weber, , Ranganath et al., ], and specifically the
implementation from
numpyro [Bradbury et al., , Bingham et al., , Phan et al., ]. We
train the variational posterior for $5000$ epochs with a learning rate of $0.01$ using the Adam optimizer. We then generate $k$ samples from the posterior.
For each algorithm, we run a sweep over all a set of hyperparameters.
We then pick the best hyperparameters and run the experiment over 10 different seeds to get the final results.
For IDS, we tune $m\in\{3,5,7,10,12\}$, $\lambda\in\{2,3,4,5\}$ and $k\in\{500,750,1000,1500\}$.
For TS and US, we tune $m\in\{3,5,7,10,12\}$ and $k\in\{500,750,1000,1500\}$.
For UCB we tune $m\in\{3,5,7,10,12\}$, $\beta\in\{0.01,0.1,0.2,0.5,1,2,5\}$ and $k\in\{500,750,1000,1500\}$, where $\beta$ controls the exploration.

Figure: Figure 8: Empirical results on the active pairwise experiments selection for the gene-gene interaction detection example using MMD-based (using RBF kernel) test statistics. The solid lines represent the mean performance, whereas the shaded region represent all the runs (min-max). IDS (ours) outperforms all the baselines significantly in terms of top scoring pairs discovered as well as the regret. In terms of known interactions, IDS still outperforms the baselines by a small margin.
Refer to caption: x14.png

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 5
presentation: 3
questions: Q1) See weaknesses section on detecting that all relevant gene pairs have been found

Q2) How many gene pairs usually have non-linear interactions? Or interactions that are significantly different from the raw pair?

Q3) Please include a baseline in which we reveal $\delta_ij$ on the pair $i,j$ for which $\delta_i, $\delta_j$ are marginally most different from $\delta_0$, perhaps taking the sum of the log density ratios.

Q4) How many samples are required for Appendix B.1 to become a faithful estimator of the densities? What happens if our densities are very far off?
soundness: 3
strengths: This paper is well-written and proposes an interesting method for active experimental design. Estimating the joint density and comparing to the marginal density is an interesting idea, both for the disjointedness analysis and for the separability analysis. 

Although prior work (e.g., GEARS) has formulated CRISPR pair perturbation response prediction as a matrix completion problem, the active learning formulation seems novel, especially in dealing with imaging modalities (as opposed to RNAseq). The authors additionally spell out their assumptions — such as a set of latent variables $Z$ which capture all the information about perturbations — clearly. Their statistical testing framework is clear and well-motivated.

The authors additionally compare to plenty of baselines on their experiments. Although they only perform experiments on one real-world dataset, experimentation in this domain is expensive, so this is a natural choice and should not be seen as a limitation of their work.

 The appendix thoroughly explains the architecture and the methodological choices.
summary: This paper presents a novel active learning method to discover pairs of highly correlated variables. Concretely, we imagine that an experimentalist has access to an observations space $\mathcal{X}$, which change as a function of perturbations. Crucially, some pairs of perturbations $\delta_i, \delta_j$ jointly operate ($\delta_{ij}$) differently on $\mathcal{X}$ relative to their marginal individual effects. The key idea of this paper is to estimate three densities — $p(x|\delta_i), p(x|\delta_j), p(x|\delta_{ij})$ — which capture the distribution of $x$ based on marginal and joint perturbations and compare them to $p(x|\delta_0)$, the world in which no perturbations are applied.

The authors formulate their active matrix completion problem in a statistical hypothesis testing framework. In particular, at each step, they estimate a posterior over $\mathbf{R}$, a matrix of test statistics for all pairs of perturbations, and greedily selects a batch of experiments to run in the next step.

The authors apply this framework to CRISPR knockouts on pairs of 50 genes. Their dataset includes observations for all pairs in this 50 gene set, so the authors simulate an active experimental design by assuming the gene pair perturbations observations are masked and then simulating matrix completion by progressively revealing observations. Their method, which they term IDS, quickly discovers top gene pairs relative to standard baselines, and minimizes natural errors like regret.
weaknesses: The use of a statistical test for matrix completion is a little poorly motivated. In particular, given an n x n matrix completion problem, at least one pair is likely to have a large test statistic. Given this is an active learning problem, multiple hypothesis testing is not as core of a concern, but I'd like the authors to discuss in their rebuttal what a "null" set of interaction pairs might look like. Concretely, suppose that in a module of 50 genes (so 1225 gene pairs) only 5 gene pairs interact with each other. Can this method be used to decide a "stopping" criterion for experimentation — i.e. would you be able to note that the posterior on the test statistics obey the null distribution once those 5 gene pairs have been found. Could you:
1. Describe the null distribution of test statistics in this framework
2. Discuss how this method could be used to decide a stopping criterion


Second, I think the authors did not explain why one would use imaging data for this matrix completion problem, instead of the more structured RNAseq data. Could you please justify why RNAseq data isn't used instead of imaging data?

Third, I did not get a sense of how often the joint distribution differs from the product of the marginals in this setting. There is some prior work that very few perturbations actually obey this relation, e.g. https://www.biorxiv.org/content/10.1101/2024.09.16.613342v1.full.pdf. How often is this method necessary, and why can't we just take experiment on the pairs $\delta_ij$ for which $\delta_i, $\delta_j$ are marginally most different from $\delta_0$. Could you provide quantitative analysis on how often joint distributions significantly differ from products of marginals in this dataset?

---

confidence: 5
contribution: 3
initial_rating: 3
presentation: 4
questions: - For research proposing a new method, it would be essential to compare its advantages and accuracy against existing approaches. There are various statistical methods (likelihood-ratio tests, two-way ANOVA, etc) to identify interaction effects, so did you conduct comparisons with these standard methods?

- In terms of detecting synthetic lethality, I believe there are several studies in genetic statistics related to identifying "epistasis." Did you consult this literature?

- What are the primary reasons you believe traditional methods are inadequate for this purpose and we need new methods?

- Two methods for testing the interaction effect are proposed, but their distinct purposes and value are unclear. Which method appears more promising, and under what conditions should each be used?

- How do you test for the KL divergence identity condition required for the "separability test"? Is this a simple binary decision based on whether the KL divergence is below a certain threshold, rather than a hypothesis test? If so, how did you determine the threshold?

- The theoretical justification for the entire Section 4 seems insufficient; could you provide additional support if available? The "interaction effect" as defined in this paper seems fundamentally based on comparing the effects of "double perturbation" with "individual perturbations." Given this, the rationale for identifying promising pairs with interaction effects using adaptive sampling or bandit without actually conducting double perturbations is unconvincing. If you view this as an application of bandit’s best-arm identification, typically, double perturbation across all (i, j) pairs would need to be measured several times to gain meaningful information on arm expectations. Could you provide the following information?
   - A clearer explanation of what prior information or assumptions allow us to predict interaction effects without measuring pairwise effects
   - More rigorous justification for the low-rank and Gaussian prior assumptions

- For the evaluation of the main claim of "automated discovery of pairwise interactions" as suggested in the paper's title, could you:
   - Clarify how the proposed approach handles the lack of repeated sampling for many (i,j) pairs?
   - Explain the assumptions and modifications needed to apply bandit methods in this setting?
   - Provide a more detailed interpretation of Figure 6, explaining what insights can be gained from the comparison to random search and conventional bandit methods?

- Have you conducted any validation beyond the synthetic lethality example? As the title suggests, this paper aims to propose a general method for testing "interaction effects," but with no examples beyond synthetic lethality, it might be more suitable to focus claims on identifying biological "epistasis."

- For instance, when generating test cases in software development, identifying bugs that occur only under certain condition combinations could fit within this paper’s defined "interaction effect." A system might function correctly when condition 1=True or condition 2=True, but malfunction when both are true. Do you think this method is also effective for such identification? Specifically, what is the rationale for identifying this condition combination as promising without observing the behavior under condition 1=True and condition 2=True? Would this not be challenging in black-box testing?
soundness: 3
strengths: - The paper addresses the problem of identifying pairwise interactions, specifically highlighting cases where the effect of two perturbations, such as cell lethality from double gene knockout, is entirely different from the effects of each perturbation alone. In the experiments, gene knockout was actually performed to validate effectiveness.
- The two proposed tests are technically intriguing. Each test is well-organized with necessary assumptions and effectively leverages existing theories, incorporating reasonable methods such as MMD-based two-sample tests (Gretton, 2012) and ASD (Xu et al., 2022) to construct a effective procedure.
- Despite potentially abstract topic of identifying "interaction effects", the Introduction clearly explains the ideas and intentions. The discussion is grounded with examples, such as validation on synthetic data and actual biological applications, making the logic easy to follow.
summary: This paper proposes a new method for exploring the structure of a black-box system of interest by introducing perturbations, specifically focusing on identifying "interacting pairs"—pairs of perturbations that yield results significantly different from the effects of each perturbation applied individually. The paper proposes two methods to examine whether given pairs of perturbations exhibit a non-trivial "interaction effect" in this sense. The first method is a "separability" test, which checks if the effects of a pair of perturbations provide information beyond that obtained from each perturbation alone. The second method is a "disjointedness" test, which quantifies whether each perturbation in the pair influences distinct subsets of the outcome space. The separability test uses a criterion based on KL divergence, while the disjointedness test uses an MMD-based two-sample test for identity. Additionally, they treat the task of calculating each test statistic for all pairs (i, j) as a matrix completion problem and apply an active learning-based sequential experimental design using ADS (Xu et al., 2022) to identify pairs (i, j) with potentially having larger test-statistics values. Experiments demonstrate the superiority of the proposed method on benchmark synthetic data that meets the study's assumptions, and they conduct a synthetic lethality test to determine if two gene knockouts result in cell lethality, empirically confirming the method’s effectiveness in real biological systems.
weaknesses: - The two proposed interaction tests are not compared with any standard methods. The problem in question is not new; it has a long history in statistics as the "interaction effect," where the combination of two or more factors produces an effect greater (or less) than the sum of their individual effects [1][2]. Traditional applied statistical methods (likelihood-ratio tests, two-way ANOVA, etc) have also been used for identifying synthetic lethality [3], so a comparative analysis and discussion of differences with conventional methods are necessary. The paper only includes internal comparisons among variations of the proposed method, offering limited basis for objectively assessing its effectiveness against previous methods.
   - [1] Cox, D. R. (1984). Interaction. International Statistical Review / Revue Internationale de Statistique, 52(1), 1–24. https://doi.org/10.2307/1403235
   - [2] Amy Berrington de González, & Cox, D. R. (2007). Interpretation of Interaction: A Review. The Annals of Applied Statistics, 1(2), 371–385. http://www.jstor.org/stable/4537441 
   - [3] Akimov Y, Aittokallio T. Re-defining synthetic lethality by phenotypic profiling for precision oncology. Cell Chem Biol. 2021 Mar 18;28(3):246-256. doi: 10.1016/j.chembiol.2021.01.026.

- Biologically, this issue is known as "epistasis," a well-researched topic with substantial existing literature. While statistical interaction tests have been applied, domain knowledge and modeling of underlying structures are essential (e.g., [4][5]). Although this study treats synthetic lethality as a special case and even conducted biological experiments, it does not discuss this background, leaving unclear what value this research or its new methods can provide on the top of many traditional studies on this long-standing research topic.
   - [4] Segrè, D., DeLuna, A., Church, G. et al. Modular epistasis in yeast metabolism. Nat Genet 37, 77–83 (2005). https://doi.org/10.1038/ng1489
   - [5] Terada A, Okada-Hatakeyama M, Tsuda K, Sese J. Statistical significance of combinatorial regulations. Proc Natl Acad Sci U S A. 2013 Aug 6;110(32):12996-3001. doi: 10.1073/pnas.1302233110.

- Consequently, the adaptive sampling and bandit approaches proposed in Section 4 to identify pairs with interaction effects “without directly measuring pairwise effects” are unconvincing and require additional justification, clearer assumptions, and validation. By definition, determining the presence of an "interaction effect" would seem to require measuring both pairwise and single effects at least. For example, in synthetic lethality, there are gene pairs where knocking out either gene alone is non-lethal, but a double knockout is lethal. In this case, we have no observable clues suggesting that this pair is likely to have an "interaction effect." While assumptions in Sections 3.1 and 3.2 seem related to this point, some, such as the "low-rank R with a Gaussian prior on the columns" in line 350, appear ad hoc and inadequately justified. This point requires thorough analysis and explanation for proper justification of the proposed method.

- Further, introducing bandit (specifically, best-arm identification) for this problem requires clarification. In bandit approaches, sampling each arm multiple times is generally assumed, making the theory inapplicable in settings with no sampling for many (i,j) pairs. For synthetic lethality, for instance, it would be necessary to test double knockouts across all (i, j) pairs with replicates. Figure 6 compares the method with random search and conventional bandit methods, but the meaning of this comparison is unclear, necessitating additional justification.

- While the paper aims to address general "interaction effect" identification, its validation is limited to synthetic lethality, providing insufficient evidence for its general effectiveness across other cases.

---

confidence: 2
contribution: 3
initial_rating: 6
presentation: 3
questions: I have the following questions.

1. By pairwise interactions do you mean causal dependencies or correlations? Can you also discover causal relations?
2. What do you mean by "unstructured data"? I would say that an image is structured as a grid of pixels.
3. Line 094 please elaborate on what you mean "by this notion of interaction".
4. Line 374 what are 3 dimensional tabular data? Please provide the dimensions.
5. Line 478 "Genes from the same pathway are ordered adjacently". Please rephrase/explain what it means.
6. Line 509. Can you explain more about how the top 5% pairs are chosen? Are they among the known interactions?
7. Line 511. How is the regret computed?

Some recommendations:
1. The first two sentences in the introduction require citation.
2. The font size in Fig. 1 is very small.
3. Increase the font size of the colormap in Fig. 2.
4. The font size in Fig. 3 (labels on images on the right) is very small.
5. In Fig. 3 you need to explain what the image on the right shows (presence or absence of interaction).
6. In Fig. 4 you should add a label in each figure and increase the font size of the colormap.
7. For Fig. 5 the same, nothing is visible.

I would like to hear your opinion on my feedback first and then be willing to raise my score.
soundness: 3
strengths: The proposed methodology is interesting. 
The main strength of this work is that it makes the biological experimental process more efficient and with lower costs. 
The theoretical claims are justified with mathematical proofs and the effectiveness of the algorithm is empirically validated.
summary: This paper proposes a method for detecting pairwise interactions between gene perturbations. 
It first proposes two statistical tests for deciding whether two perturbations are separable or disjoint.
The statistical test of disjointness is employed in a greedy active matrix completion algorithm that decides which pair of perturbations to examine in the next step. The proposed method is evaluated in sythetical settings showing the validity of the statistical tests as well as in a real gene perturbation experiment where it performs better against other baselines that consider different choice of policy in the active matrix completion procedure.
weaknesses: I am not a specialist in the biological field of gene perturbations experiments, but based on my understanding I would point out the following potential weaknesses for the improvement of the paper.

1. It is not quantified how much the biological experiments benefit from the active learning algorithm in terms of the total number of necessary perturbations.
How much more efficient is your algorithm compared to a standard exhaustive approach that would consider all possible perturbations? 

2. I am not sure if it possible but it would make sense to additionally compare the proposed method against causal discovery methods that utilize interventions, for example, DCDI [1] or GIES [2], or methods specifically designed for gene regulatory networks [3]. It would be interesting to see what would be the performance of a causal discovery method for a subset of given perturbations. Ofcourse that would require transforming the data into a format that contains a single value for every node for each sample.

[1] Brouillard, Philippe, et al. "Differentiable causal discovery from interventional data." Advances in Neural Information Processing Systems 33 (2020): 21865-21877.

[2] Hauser, Alain, and Peter Bühlmann. "Characterization and greedy learning of interventional Markov equivalence classes of directed acyclic graphs." The Journal of Machine Learning Research 13.1 (2012): 2409-2464.

[3] Aibar, Sara, et al. "SCENIC: single-cell regulatory network inference and clustering." Nature methods 14.11 (2017): 1083-1086.

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 3
questions: * Can you counter my points (in weaknesses) regarding the lack of novelty in this work?

Minor points:
* The right-hand side of the bottom equation on p.4 seems off; isn't $p_{Z_i}$ a distribution over a scalar random variable, and the argument $g^{-1}(x)$ is a vector?
* I think the clarity can be improved (and some space saved) if you apply a log to both sides of Equation (1) and go directly to Equation (3), without writing Equation (2).
soundness: 3
strengths: * The authors tackle the important problem of determining whether it is worthwhile to run a costly experiment where two perturbations are applied simultaneously.
* The evaluation is thorough, and demonstrates utility on a real-world problem.
summary: The authors propose two test statistics to measure the separability and disjointness between pairs of perturbations. Perturbations are separable if they affect disjoint subsets of latent variables, and are disjoint if they affect disjoint parts of the observation. Their separability metric is validated on two synthetic datasets where the ground-truth data generating process is known (i.e. the graph of the perturbation and latent variables), as well as on a real cell painting dataset. They also show that the disjointness metric can be used in an active learning context on cell painting data to better uncover perturbations with pairwise interactions relative to several baselines.
weaknesses: * This is an applied paper, and I believe the novelty is weak.
* The separability metric essentially measures the mutual information between the latent variables corresponding to a pair of perturbations. There are many different ways to instantiate this, and it's unclear whether the authors' proposed metric itself is a valuable contribution.
* A similar statement also holds for the disjointness metric.
* Essentially, this paper combines many existing methods to achieve an empirical (and somewhat narrow) goal. I believe this is a valuable contribution, but it may be suited for a more applied venue.

Meta-review: An interesting paper, which unfortunately does not pass the bar for acceptance at ICLR. Substantial points of disagreements still coexisted with 1aqQ at the end of the review process, and while I can buy the authors' argument for point (1), I believe a bit more could have been done for points (2-3). I can only recommend that the authors polish further their paper for a next round of submission.

Decision: Reject