#  RECOVER: sequential model optimization platform for combination drug
repurposing identifies novel synergistic compounds in vitro

Paul Bertin1, Jarrid Rector-Brooks1, Deepak Sharma1, Thomas Gaudelet2, Andrew
Anighoro2, Torsten Gross2, Francisco Martínez-Peña3, Eileen L. Tang3, Suraj M
S2, Cristian Regep2, Jeremy B. R. Hayter2, Maksym Korablyov1, Nicholas
Valiante4, Almer van der Sloot5, Mike Tyers5, Charles Roberts2, Michael M.
Bronstein6,7, Luke L. Lairson3, Jake P. Taylor-King2, and Yoshua Bengio1
1 Mila, the Quebec AI Institute, Canada 2 Relation Therapeutics, London, UK
3 Department of Chemistry, The Scripps Research Institute, USA 4 Glyde Bio,
Inc, USA.
5 IRIC, Institute for Research in Immunology and Cancer, Université de
Montréal, Canada
6 Department of Computer Science, University of Oxford, UK 7 Twitter, UK
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

For large libraries of small molecules, exhaustive combinatorial chemical
screens become infeasible to perform when considering a range of disease
models, assay conditions, and dose ranges. Deep learning models have achieved
state of the art results in silico for the prediction of synergy scores.
However, databases of drug combinations are biased towards synergistic agents
and these results do not necessarily generalise out of distribution. We employ
a sequential model optimization search utilising a deep learning model to
quickly discover synergistic drug combinations active against a cancer cell
line, requiring substantially less screening than an exhaustive evaluation.
After only $3$ rounds of ML-guided in vitro experimentation (including a
calibration round), we find that the set of drug pairs queried is enriched for
highly synergistic combinations; two additional rounds of ML-guided
experiments were performed to ensure reproducibility of trends. Remarkably, we
rediscover drug combinations later confirmed to be under study within clinical
trials. Moreover, we find that drug embeddings generated using only structural
information begin to reflect mechanisms of action.

Our method is available at: https://github.com/RECOVERcoalition/Recover.

##  1\. Introduction

Figure 1. Overview of the RECOVER workflow integrating both a novel machine
learning pipeline and iterated wet lab evaluation.

Drug combinations are an important therapeutic strategy for treating diseases
that are subject to evolutionary dynamics, in particular cancers and
infectious disease (tyers2019drug, ; mokhtari2017combination, ). Conceptually,
as tumours or pathogens are subject to change over time, they may develop
resistance to a single agent (delou2019highlights, ) — motivating one to
target multiple biological mechanisms simultaneously (al-
lazikaniCombinatorialDrugTherapy2012, ). Discovering synergistic drug
combinations is a key step towards developing robust therapies as they hold
the potential for greater efficacy whilst reducing dose, and hopefully
limiting the likelihood of adverse effects. For example, in a drug repurposing
scenario (i.e., uncovering new indications for known drugs), the ReFRAME
library of $\sim$12,000 clinical stage compounds (janes2018reframe, ) leads to
$\sim$72 million pairwise combinations; this does not appear tractable with
standard high throughput screening (HTS) technology — even at a single dose
(clare2019industrial, ).

With the recent COVID-19 global health crisis, there has been the need for
rapid drug repurposing that would allow for expedited and de-risked clinical
trials. Due to the complexity of selecting drug combinations and minimal
training data publicly available, studies have typically been limited towards
monotherapy repurposing from a variety of angles — often involving artificial
intelligence (AI) techniques to provide recommendations (zhou2020artificial,
). The dearth of drug combination datasets is due to the large combinatorial
space of possible experiments available — ultimately limiting the quality of
drug synergy predictions. Moreover, databases of drug combinations are biased
towards suspected synergistic agents and thus making predictions outside the
scope of the training dataset can be challenging. To avoid requiring large
amounts of data, an emerging cost-efficient solution involves the use of
sequential model optimization (SMO). Within SMO, one performs both highly
informative experiments (“exploration”) and experiments that double-down on a
promising hypothesis (“exploitation”) (sverchkov2017review, ). Utilisation of
such approaches can result in discovering most of the promising combinations
within a given library whilst requiring substantially less experimentation
than an exhaustive search. With an SMO platform available in conjunction with
an appropriate in vitro assay, one has a powerful tool to rapidly respond to a
future public health crisis.

There have now been a number of approaches for predicting the effects of and
subsequently prioritising drug combinations
(bulusuModellingCompoundCombination2016, ). Classic bioinformatics approaches
have focused on using machine learning and network statistics over specified
features of drugs (e.g., molecular fingerprints (cereto2015molecular, )), cell
lines (e.g., transcriptomics, copy number variations
(karczewski2018integrative, )), and interactome topology between biomolecules
(e.g., protein–protein interactions, chemical–genetic interactions
(wildenhain2015prediction, ), or gene regulatory networks (cheng2019network,
)). Initiatives such as the Dialogue on Reverse Engineering Assessment and
Methods (DREAM) have led to a plethora of methods being benchmarked against
one another in prospective challenges through the generation of novel datasets
(menden2019community, ). Complex deep learning architectures, which have set
state of the art performance across a number of domains (lecun2015deep, ),
have been used to predict both adverse drug–drug interactions
(zitnik2018modeling, ; deac2019drug, ) and synergistic drug combinations
(preuer2018deepsynergy, ; jin2021deep, ; rozemberczki2021moomin, ). Sequential
approaches, wherein several rounds of selection are performed, have also been
explored in the context of drug combinations; for example, Kashif et al.
(Kashif_Andersson_Hassan_Karlsson_Senkowski_Fryknas_Nygren_Larsson_Gustafsson_2015,
) have proposed a heuristic based (as opposed to model based) exploration
strategy.

We present a sequential model optimization platform that can guide wet lab
experiments: RECOVER. We first propose a deep learning regression model that
predicts synergy using molecular fingerprints as inputs. We perform a
retrospective validation to benchmark the performance of our model and
understand its generalization abilities using the DrugComb database — largely
pertaining to cancer cell line data (zagidullin2019drugcomb, ). Thereafter, we
evaluate our SMO pipeline in silico, which allows the model to select the most
relevant data points to be labelled in order to discover most promising
combinations while reducing model uncertainty. Finally, we test RECOVER
prospectively in an in vitro experimental setting whereby we discover novel
synergistic combinations active against a breast cancer model cell line, MCF7,
which is well represented within our training dataset. This work forms a proof
of concept demonstration of RECOVER — which should then motivate greater
community adoption of the method and extensions thereof.

##  2\. Results

###  2.1. RECOVER: sequential model optimization platform for rapid drug
repurposing

RECOVER is an open-source SMO platform for the optimal suggestion of drug
combinations, see Figure 1. To encourage use by the scientific community, we
detail a configuration that can be trained on a personal computer or laptop
without requiring dedicated computational infrastructure. Pairs of drug
feature vectors are fed into a deep neural network which is used for the
prediction of synergy scores. These feature vectors include molecular
fingerprints as well as a one-hot encoding identifying a drug. For a
description of the model, see Appendix A.1.

Our core focus is the prediction of pairwise drug combination synergy scores.
Whilst many mathematical descriptions of synergy have been proposed
(tyers2019drug, ), in the following work, we utilise the Bliss synergy score
due to its simplicity and numerical stability. In the context of cell
viability, the Bliss independence model assumes that in the absence of
synergistic effects, the expected fraction of viable cells after treatment
with drugs $a$ and $b$ at doses $c_{a}$ and $c_{b}$, written $V(c_{a},c_{b})$,
is identical to the product of the fractions of viable cells when utilising
each drug independently, i.e., $V(c_{a})V(c_{b})$. We then define the Bliss
synergy score as the difference between these quantities such that a fraction
of surviving cells $V(c_{a},c_{b})$ smaller than the expected proportion
$V(c_{a})V(c_{b})$ leads to a large Bliss synergy score

(1) $\displaystyle\textnormal{s}_{Bliss}(c_{a},c_{b})$
$\displaystyle=V(c_{a})V(c_{b})-V(c_{a},c_{b})$
$\displaystyle=I(c_{a},c_{b})-I(c_{a})-I(c_{b})+I(c_{a})I(c_{b}),$

where $I(\cdot)=1-V(\cdot)$ is the experimentally measured growth inhibition
induced by drug $a$, $b$, or both together at the associated doses. Given a
dose-response matrix for the two drugs, a global synergy score can be obtained
through a pooling strategy. In our case, we take the maximum value, i.e.,

(2)
$\displaystyle\hat{\textnormal{s}}_{Bliss}=\max_{c_{a},c_{b}}\,\textnormal{s}_{Bliss}(c_{a},c_{b})\,.$

In many studies, the arithmetic mean is taken to calculate a global synergy
score. Unfortunately, different laboratories use different dose intervals for
each drug, and typically each drug combination shows a synergistic effect at a
specific dose-pair interval. Therefore, the arithmetic mean is highly
sensitive to the chosen dose interval, thus why we choose to prioritise a max-
pooling strategy as in Eq. 2. Unless explicitly stated otherwise, a synergy
score refers to a global max-pooled Bliss score.

In addition to the prediction of synergy, RECOVER estimates the uncertainty
associated with the underlying prediction. More precisely, for a given
combination of drugs, RECOVER does not only provide a point estimate of the
synergy, but estimates the distribution of possible synergy scores for each
combination, which we refer to as the predictive distribution. We define the
uncertainty as the standard deviation of the predictive distribution.

An acquisition function is used to select the combinations that should be
tested in subsequent experiments (vzilinskas1972bayes, ). This acquisition
function is designed to balance between: exploration, prioritising
combinations with high uncertainty whereby labelling said points should
increase predictive accuracy in future experimental rounds; and exploitation,
selection of combinations believed to be synergistic with high confidence.

In summary, this SMO setting consists of generating recommendations of drug
combinations that will be tested in vitro at regular intervals. At each step,
RECOVER is trained on all the data acquired up to that point, and predictions
are made for all combinations that could be hypothetically tested
experimentally. The acquisition function is then used to provide
recommendations for in vitro testing. The results of the experiments are then
added to the training data for the next round of experiments and the whole
process repeats itself.

###  2.2. Retrospective testing of RECOVER informs the design of future
experiments

We evaluate the in silico performance of RECOVER in 3 ways. We first compare
against baselines, and evaluate various aspects of the performance of RECOVER
(see Appendix B). Then, and most crucially, we test its ability to generalise
in several out-of-distribution tasks in order to understand the scope of
scenarios to which it can be applied. This analysis is the focus within this
section. Finally, we design our experimental setting based on these insights
and perform backtesting through simulating mock SMO experiments (see
Appendices A.2 and C).

Due to the limited size of most individual drug combination studies reported
in the literature, we focus on the NCI-ALMANAC viability screen
(holbeck2017national, ) summarised in Figure S2. Unless explicitly stated
otherwise, presented in silico benchmarking utilises this resource. We refrain
from combining multiple datasets due to the severe batch effects between
studies; in Figure S2E, we show a scatter plot that demonstrates inconsistency
between the O’Neil 2016 (o2016unbiased, ) series of drug combination
experiments against their NCI-ALMANAC counterpart. We note this may result
from: variation in the readouts of these experiments, mutations in cell lines,
or differences in harvest times.

We investigate whether RECOVER can generalize beyond the training distribution
in various ways: (a.) What is the performance on test cases drawn from the
same distribution as the training set? (b.) Can RECOVER generalize to unseen
drug pairs? (c.) Can RECOVER leverage observations from other cell lines in
order to improve the predictions on the MCF7 cell line? (d.) Can RECOVER
predict the synergy of an unseen combination on the MCF7 cell line (using
information on the drug combination tested in other cell lines)? (e.) Can
RECOVER predictions transfer between the NCI-ALMANAC and O’Neil 2016 studies?
These tasks are illustrated graphically in Figure 2. Through understanding the
capacity for the generalisation of RECOVER, we can design prospective
experiments with greater confidence of success.

Depending on the task at hand, the model configuration will differ slightly.
For all tasks (excluding task (b.)), the drug feature representations consists
of the Morgan fingerprints (morgan1965generation, ) concatenated with a one-
hot encoding specifying the identity of drug. For task (b.), only the Morgan
fingerprint is used. When several cell lines are available, the model is
conditioned on cell lines using feature-wise linear modulation (FiLM)
(perez2018film, ), see Appendix A.1.3. The cell line features are either a
one-hot encoding of the cell line for tasks (c.) and (e.), or some information
about mutations and basal level of mRNA gene expression for task (d.).

In Table 1, we report the validation and test performance metrics of RECOVER
across each of the five tasks. Examining performance within task (a.), test
statistics appear modest, however we demonstrate limits on achievable
performance in Appendix B.3 — resulting from experimental noise and non-
uniformity of synergy scores. In task (b.) we note a considerable drop in
performance when compared to task (a.), demonstrating RECOVER will have
markedly reduced performance when attempting to predict the synergy of drug
combinations when both drugs have not been observed in earlier experiments.
The results pertaining to tasks (c.) and (d.) demonstrate that leveraging
experiments from other cell lines does provide a benefit when compared to the
performance from task (a.), although the effect is most significant when the
specific drug combination in question has been seen in other cell lines, i.e.,
task (d.). For completion, we confirm the significant batch effects between
the NCI-ALMANAC and the O’Neil 2016 studies render using the same model
parameters for both studies impossible — notice task (e.) performing at the
level of randomness and Figure S2E.

From the above results, we can recommend that any prospective experiments
should require that one of the two drugs in the combination have been seen in
some context before, see task (b.). In addition, should we want to utilise
publicly available resources, we will have to intelligently incorporate
existing publicly available resources, see task (e.). To this end we
investigated using transfer learning, wherein one trains a model on a large
dataset (known as pretraining), and thereafter refines the model on a smaller
dataset (known as fine-tuning) — typically with some aspect of the task or the
data changed between the two instances. We show that this is possible and
beneficial (compared to not leveraging existing data) in an SMO setting
between the O’Neil 2016 and NCI-ALMANAC studies in Appendix C.3. Remarkably
even with minimal correlation between studies, we are able to observe the
benefits of transfer learning in this scenario. We rely upon these findings to
make design choices regarding our prospective experimental work.

\begin{overpic}[width=433.62pt]{img/main/Figure2.pdf} \put(3.0,-1.5){(a.) Default} \put(15.0,-1.5){(b.) Drug level split} \put(34.0,-1.5){(c.) Multi cell line} \put(53.0,-1.5){(d.) Cell line transfer} \put(76.0,-1.5){(e.) Study transfer} \end{overpic} Figure 2. Overview of the different tasks on which RECOVER has been evaluated. Each task corresponds to a different way to split the training, validation and test sets, and aims at evaluating a specific generalization ability of the model. (a.) Default. Combinations are split randomly into training/validation/test (70%/20%/10%). Only the MCF7 cell line is used. (b.) Drug Level Split. The test set consists of combinations among drugs that are excluded from the training and validation sets. 30% of available drugs are used for test. Remaining combinations are split into training and validation (80%/20%). Only the MCF7 cell line is used. (c.) Multi Cell Line. All cell lines are used for training and validation, but the test set is restricted to the MCF7 cell line. If a combination is part of the test set, it is never used in the training and validation sets, regardless of cell line. Training/validation split is consistent across cell lines. Proportions (70%/20%/10%) (d.) Cell Line Transfer. The test set consists of all examples corresponding to the MCF7 cell line. Other cell lines are randomly assigned to training/validation (80%/20%). A given combination can appear in both training and test sets (e.) Study Transfer. The O’Neil 2016 study (o2016unbiased, ) is used to generate the training and validation sets. All overlapping cell lines are used. The test set is generated using NCI-ALMANAC and restricted to combinations for which both drugs also appear in the O’Neil study. If a combination is part of the test set, it is excluded from training and validation sets, regardless of cell line. |  | (a.) Default | (b.) Drug Level Split | (c.) Multi Cell Line | (d.) Cell Line Transfer | (e.) Study Transfer
---|---|---|---|---|---|---
$R^{2}$ | Valid. | $0.343\pm 0.053$ | $0.401\pm 0.147$ | $0.387\pm 0.032$ | $0.278\pm 0.027$ | $0.304\pm 0.021$
Test | $0.242\pm 0.006$ | $0.038\pm 0.002$ | $0.282\pm 0.017$ | $0.382\pm 0.017$ | $0.014\pm 0.016$
Spearman c. | Valid. | $0.474\pm 0.021$ | $0.459\pm 0.069$ | $0.518\pm 0.021$ | $0.299\pm 0.047$ | $0.589\pm 0.032$
Test | $0.466\pm 0.007$ | $0.157\pm 0.012$ | $0.448\pm 0.021$ | $0.378\pm 0.015$ | $0.147\pm 0.075$
Table 1. Performance of RECOVER for the different tasks, as detailed in Figure
2. Standard deviation computed over 3 seeds.

###  2.3. Prospective use of RECOVER enriches for selection of synergistic
drug combinations

\begin{overpic}[trim=2.84544pt 2.84544pt 2.84544pt
2.84544pt,clip,width=182.1196pt]{img/main/Figure3A.png} \put(-5.0,80.0){\Large
A} \put(30.0,70.0){\footnotesize 149 single drugs}
\put(30.0,63.0){\footnotesize 2,916 available combinations} \end{overpic}

\begin{overpic}[clip,width=242.82613pt]{img/main/Figure3B_ECDF.png}
\put(52.75,29.0){\includegraphics[clip,width=106.69783pt]{img/main/Figure3B_Boxplot.png}}
\put(38.0,12.0){\includegraphics[clip,width=78.24507pt]{img/main/Figure3B_Calibration_Round.png}}
\put(54.0,27.5){\Huge{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{/}}}
\put(1.0,60.0){\Large B} \end{overpic}

Figure 3. Prospective use of RECOVER for in vitro evaluation. (A.) Heatmap
representing drug combinations used during pretraining (NCI-ALMANAC), in five
subsequent experiment rounds, and excluded from the analysis. Drug
combinations which were not available for pretraining or were not selected for
experiments are represented in white. (B.) Cumulative density plot of max
Bliss synergy score for each experimental round; (inset) box plot
representation and calibration round details.

From the in silico results, we now test RECOVER prospectively using a cancer
cell model, leveraging publicly available data for pretraining. Using the
insights from Section 2.2, the queriable space of drug combinations was
designed to include drug pairs where only one compound was already seen by the
model during pretraining — with a second compound not seen before. Details
about the model used to generate recommendations are available in Appendix
A.3.

The MCF7 cell line was used to generate $6\times 6$ dose-response matrices,
see Appendix A.5 for more details. We perform multiple rounds of RECOVER-
informed wet lab experiments and observe sequential improvements in
performance. The rounds of experiments are described as follows:

  * 1\. Calibration.

The initial round of experiments was performed to supplement publicly
available data with 20 randomly selected unseen drug combinations.
Furthermore, we confirmed the previous in silico result that we could not
predict synergy scores for unseen drugs (using public data) through selecting
5 highly synergistic combinations selected by RECOVER. In addition, 5 more
drug combinations were selected by a Graph Neural Network (GNN) model in the
style of Zitnik et al. (zitnik2018modeling, ) that we did not develop further
due to the computational overhead. It was also specified that each drug should
appear in at most a single drug combination queried.

  * 2\. Diversity.

Thereafter, drug combinations are selected using using model predictions in
conjunction with the upper confidence bound (UCB) acquisition function. To
ensure that we quickly observe all single drugs at least once (as we showed
that the model cannot generalize well to combinations involving unseen drugs),
we select our batch of experiments as follows. First we rank combinations
according to their acquisition function score. We then find the first
combination which involves a drug that has not yet been used (or involved in
one of the combinations from the current batch), and add it to the batch. We
repeat until we have 30 combinations in the batch.

  * 3\. SMO Search.

RECOVER is now free to select any drug pairs of interest for testing, with the
requirement that any single drug may be selected no more than 5 times (to
avoid oversampling and depletion of chemical stock). Three such rounds have
been performed in this manner.

The search space was constructed as follows. The NCI-ALMANAC includes 95
unique drugs that were employed in combinations tested on the MCF7 cell line,
see grey area in Figure 3A. We chose to deprioritise drugs without a well-
characterised mechanism of action (MoA) to facilitate biological
interpretation and validation of the results, see light blue area in Figure
3A. To achieve this, drugs in NCI-ALMANAC were annotated with known targets
extracted from the ChEMBL drug mechanism table: 54 drugs matched with at least
one known target were thus selected. An additional 54 drugs were selected by
clustering drugs with known MoA that are included in the DrugComb
(zagidullin2019drugcomb, ) database but not in NCI-ALMANAC, see Appendix A.4
for details. Hence, a search space including a total of 2,916 drug
combinations was obtained, see the white area in the left bottom corner of
Figure 3A.

In Figure 3B, we plot the cumulative density function of each experimental
round. We note that the mean synergy score significantly increases between the
first and the third round (t-test, $\text{p}<0.05$); this trend further
continues by the fifth round (t-test, $\text{p}<10^{-5}$). Moreover, the
distribution starts developing a heavier tail towards high max Bliss synergy
scores. This emergent heavy tail already appears significant when comparing
the distribution in the first SMO Search round to the background distribution
of synergy scores in NCI-ALMANAC (Kolmogorov–Smirnov test, $\text{p}<0.025$),
see Figure S2B. Finally, the highest max Bliss synergy score observed
increases between rounds until the second SMO Search round, whereby the
behaviour appears to have saturated.

###  2.4. Discovery and rediscovery of novel synergistic drug combinations

In Figure S3, we plot the Bliss synergy dose-response matrices for the 12 most
synergistic drug combinations (from $\sim$150 tested) with Alfacalcidol and
Crizotinib achieving a max Bliss score above 90. Of note, we rapidly discover
drug combinations with similar mechanisms and efficacy to those already in
clinical trials. Namely, within the first SMO search round we found: (a.)
Alisertib & Pazopanib, and (b.) Flumatinib & Mitoxantrone. Compared to the
former combination, in the latter combination neither drug requires a high
dose of 50 $\mu$M to achieve the maximum Bliss synergy score.

We plot the Bliss synergy matrices in Figure S4. Pazopanib inhibits
angiogenesis through targeting a range of kinases including vascular
endothelial growth factor receptor (VEGFR), platelet-derived growth factor
receptor (PDGFR), c-KIT and fibroblast growth factor receptors (FGFR); in
contrast Alisertib is a highly selective inhibitor of mitotic Aurora A kinase.
Synergism between the two agents is hypothesised to be linked to the
observation that mitosis-targeting agents also demonstrate antiangiogenic
effects. In an independent study, the combination of Alisertib & Pazopanib has
successfully completed phase 1b clinical trials for advanced solid tumours
(shah2019phase, ). The combination of Flumatinib and Mitoxantrone appears to
be linked to a similar mechanism, but does not seem to have been studied in
the biomedical literature. Whilst Flumatinib is a tyrosine kinase inhibitor
targeting Bcr-Abl, PDGFR and c-KIT, Mitoxantrone is a Type II topoisomerase
inhibitor.

###  2.5. RECOVER drug embeddings capture both structural and biological
information

\begin{overpic}[trim=0.28436pt 0.28436pt 0.28436pt
0.28436pt,clip,width=429.28616pt]{img/main/Figure4.png} \end{overpic} Figure
4. UMAP of RECOVER drug embeddings with the colour scheme generated to
indicate the known target profile of the drugs. .

To get a better insight into the drug embeddings learnt by RECOVER, we report
uniform manifold approximation and projection (UMAP) visualizations of the
drug embeddings generated by the single drug module in Figure 4. The colour of
each point is chosen by applying principal component analysis (PCA) to the
binary matrix of drug–targets and scaling the first 3 dimensions into an RGB
triplet; high transparency indicates drugs with a PCA target profile close to
the average PCA target profile (calculated over all drugs).

We note that the RECOVER model does not use information on drug targets,
however drugs with similar colours are located within similar areas of UMAP
space. We also observe broad sensible patterns in UMAP space based on
structure, for example, most kinase inhibitors (with the -nib suffix) appear
in the upper left hand of the UMAP. Moreover, drugs with similar mechanisms
tend to be co-located, for example, see structurally diverse DNA-targeting
agents in the bottom right of the UMAP. As a counterpoint, we observe that
agents with either: mixed agonist/antagonist profiles, including selective
estrogen receptor modulators (SERMS); or targeting genes through indirect
mechanisms, including mammalian target of rapamycin (mTOR), lead to less
structured patterns in UMAP space. We believe this is a highly novel
observation, and suggests that were this screen to be scaled to a larger
library of small molecules, one may be able to group diverse structures into
common biological mechanisms.

##  3\. Discussion

Drug combinations can achieve benefits unattainable by mono-therapies and are
routinely investigated within clinical trials (e.g., PD-1/PD-L1 inhibitors
combined with other agents (zhu2021combination, )), and utilised within
clinical practice (e.g., antiretroviral treatment of HIV where between 3-4
agents may be used (feng2019quadruple, )). To this end, we have presented the
SMO toolbox RECOVER for drug combination identification. We showcase a general
methodology, consisting of careful analysis of the properties of our machine
learning pipeline — such as its out-of-distribution generalization capacities
— to help us design key aspects of our prospective experiments, to eventually
ensure a smooth and successful interaction between the SMO pipeline and the
wet lab. Highly synergistic drug combinations have been identified and the
resulting learnt embeddings appear to capture both structural and biological
information. We provide commentary on key aspects on our approach covering
datasets, computational methodology, and wet-lab techniques.

We note the considerable difficulties of working with publicly available
datasets with discrepancies in the data generation process. Inconsistent media
between multiple labs, the presence of spontaneous mutations within
immortalised in vitro models, and differences in experimental protocols limit
ease of data integration between laboratories (hirsch2019vitro, ). In
particular, systematic biases limit generalisability of model predictions to
subsequent prospective experiments. Within oncology, protein coding mutations
may drive resistance to any one chemotherapeutic agent, but also large scale
gene dosing changes from non-coding mutations (avsec2021effective, ), copy
number variation (shao2019copy, ) and aneuploidy (taylor2021rethinking, ).
These issues have been somewhat alleviated through careful choice of metric to
optimise (max pooled Bliss synergy scores) and only using publicly available
data for pretraining (when compared to using this data for prediction without
adaptation).

From a computational perspective, we experimented with a range of more
complicated models. For example, we considered using graph neural networks to
model biomolecular interactions (gaudelet2021utilizing, ), which have numerous
benefits including greater biological interpretability and incorporation of
prior knowledge, namely drug–target and protein–protein interactions. However,
these models only resulted in marginal increases in performance whilst
requiring substantially more computational resources. We believe that the
limited diversity of the dataset and the simplicity of the task, a one
dimensional regression, did not allow these more advanced approaches to reach
their full potential. Therefore, we prioritised a framework that could be run
quickly for rapid turnaround of recommendations for experimental testing.

When considering a SMO setting, we are required to collapse highly complex
information into a single number to be optimised (i.e., a synergy score).
Whilst there is opportunity to improve choices of metric (synergy scores may
not reflect absolute cell viability), assay readouts that better characterise
cell state (compared to cell viability) may provide a stronger starting point.
In particular, an omic readout, through transcriptomics (subramanian2017next,
) and/or single cell profiling (chen2021high, ; peidli2022scperturb, ); and
high content imaging (bray2016cell, ) provide a much higher dimensional
measurement of cell state. Furthermore, derived properties from these readouts
may be more interpretable, e.g., pathway activation (nguyen2019identifying, )
or extracellular signalling (taylor2019simulated, ). Remarkably, even while
only using cell viability as a readout, we achieved significant progress in
identifying novel synergistic drug combinations.

From the systematic screen by Jaak et al. (jaaks2022effective, ), they
conclude that: synergy between drugs is rare and highly context-dependent.
RECOVER provides a means to identify such synergies while requiring
substantially less screening than an exhaustive evaluation; thus we expect
RECOVER and similar such systems may have a role to play when addressing novel
emergent infectious disease such as the COVID-19 pandemic.

##  4\. Acknowledgements

This work was supported, in whole or in part, by the Bill & Melinda Gates
Foundation [INV-019229]. Under the grant conditions of the Foundation, a
Creative Commons Attribution 4.0 Generic License has already been assigned to
the Author Accepted Manuscript version that might arise from this submission.
The findings and conclusions contained within are those of the authors and do
not necessarily reflect positions or policies of the Bill & Melinda Gates
Foundation. The authors would also like to thank Andrew Trister and Isabelle
Lacroix for useful discussion and support.

## References

  * [1] Mike Tyers and Gerard D Wright.  Drug combinations: a strategy to extend the life of antibiotics in the 21st century.  Nature Reviews Microbiology, 17(3):141–155, 2019.
  * [2] Reza Bayat Mokhtari, Tina S Homayouni, Narges Baluch, Evgeniya Morgatskaya, Sushil Kumar, Bikul Das, and Herman Yeger.  Combination therapy in combating cancer.  Oncotarget, 8(23):38022, 2017.
  * [3] João Delou, Alana SO Souza, Leonel Souza, and Helena L Borges.  Highlights in resistance mechanism pathways for combination therapy.  Cells, 8(9):1013, 2019.
  * [4] Bissan Al-Lazikani, Udai Banerji, and Paul Workman.  Combinatorial drug therapy for cancer in the post-genomic era.  Nature Biotechnology, 30(7):679–692, July 2012.
  * [5] Jeff Janes, Megan E Young, Emily Chen, Nicole H Rogers, Sebastian Burgstaller-Muehlbacher, Laura D Hughes, Melissa S Love, Mitchell V Hull, Kelli L Kuhen, Ashley K Woods, et al.  The reframe library as a comprehensive drug repurposing library and its application to the treatment of cryptosporidiosis.  Proceedings of the National Academy of Sciences, 115(42):10750–10755, 2018.
  * [6] Rachel H Clare, Catherine Bardelle, Paul Harper, W David Hong, Ulf Börjesson, Kelly L Johnston, Matthew Collier, Laura Myhill, Andrew Cassidy, Darren Plant, et al.  Industrial scale high-throughput screening delivers multiple fast acting macrofilaricides.  Nature communications, 10(1):1–8, 2019.
  * [7] Yadi Zhou, Fei Wang, Jian Tang, Ruth Nussinov, and Feixiong Cheng.  Artificial intelligence in covid-19 drug repurposing.  The Lancet Digital Health, 2020.
  * [8] Yuriy Sverchkov and Mark Craven.  A review of active learning approaches to experimental design for uncovering biological networks.  PLoS computational biology, 13(6):e1005466, 2017.
  * [9] Krishna C. Bulusu, Rajarshi Guha, Daniel J. Mason, Richard P. I. Lewis, Eugene Muratov, Yasaman Kalantar Motamedi, Murat Cokol, and Andreas Bender.  Modelling of compound combination effects and applications to efficacy and toxicity: State-of-the-art, challenges and perspectives.  Drug Discovery Today, 21(2):225–238, February 2016.
  * [10] Adrià Cereto-Massagué, María José Ojeda, Cristina Valls, Miquel Mulero, Santiago Garcia-Vallvé, and Gerard Pujadas.  Molecular fingerprint similarity search in virtual screening.  Methods, 71:58–63, 2015.
  * [11] Konrad J Karczewski and Michael P Snyder.  Integrative omics for health and disease.  Nature Reviews Genetics, 19(5):299, 2018.
  * [12] Jan Wildenhain, Michaela Spitzer, Sonam Dolma, Nick Jarvik, Rachel White, Marcia Roy, Emma Griffiths, David S Bellows, Gerard D Wright, and Mike Tyers.  Prediction of synergism from chemical-genetic interactions by machine learning.  Cell Systems, 1(6):383–395, 2015.
  * [13] Feixiong Cheng, István A Kovács, and Albert-László Barabási.  Network-based prediction of drug combinations.  Nature communications, 10(1):1–11, 2019.
  * [14] Michael P Menden, Dennis Wang, Mike J Mason, Bence Szalai, Krishna C Bulusu, Yuanfang Guan, Thomas Yu, Jaewoo Kang, Minji Jeon, Russ Wolfinger, et al.  Community assessment to advance computational prediction of cancer drug combinations in a pharmacogenomic screen.  Nature communications, 10(1):1–17, 2019.
  * [15] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [16] Marinka Zitnik, Monica Agrawal, and Jure Leskovec.  Modeling polypharmacy side effects with graph convolutional networks.  Bioinformatics, 34(13):i457–i466, 2018.
  * [17] Andreea Deac, Yu-Hsiang Huang, Petar Veličković, Pietro Liò, and Jian Tang.  Drug-drug adverse effect prediction with graph co-attention.  arXiv preprint arXiv:1905.00534, 2019.
  * [18] Kristina Preuer, Richard PI Lewis, Sepp Hochreiter, Andreas Bender, Krishna C Bulusu, and Günter Klambauer.  Deepsynergy: predicting anti-cancer drug synergy with deep learning.  Bioinformatics, 34(9):1538–1546, 2018.
  * [19] Wengong Jin, Jonathan M Stokes, Richard T Eastman, Zina Itkin, Alexey V Zakharov, James J Collins, Tommi S Jaakkola, and Regina Barzilay.  Deep learning identifies synergistic drug combinations for treating covid-19.  Proceedings of the National Academy of Sciences, 118(39), 2021.
  * [20] Benedek Rozemberczki, Anna Gogleva, Sebastian Nilsson, Gavin Edwards, Andriy Nikolov, and Eliseo Papa.  Moomin: Deep molecular omics network for anti-cancer drug combination therapy.  arXiv preprint arXiv:2110.15087, 2021.
  * [21] M. Kashif, C. Andersson, S. Hassan, H. Karlsson, W. Senkowski, M. Fryknäs, P. Nygren, R. Larsson, and M.G. Gustafsson.  In vitro discovery of promising anti-cancer drug combinations using iterative maximisation of a therapeutic index.  Scientific Reports, 5(1):14118, Nov 2015.
  * [22] Bulat Zagidullin, Jehad Aldahdooh, Shuyu Zheng, Wenyu Wang, Yinyin Wang, Joseph Saad, Alina Malyutina, Mohieddin Jafari, Ziaurrehman Tanoli, Alberto Pessia, et al.  Drugcomb: an integrative cancer drug combination data portal.  Nucleic acids research, 47(W1):W43–W51, 2019.
  * [23] A Žilinskas and J Mockus.  On a bayes method for seeking an extremum.  Automatika i vychislitelnaja tekhnika, (3), 1972.
  * [24] Susan L Holbeck, Richard Camalier, James A Crowell, Jeevan Prasaad Govindharajulu, Melinda Hollingshead, Lawrence W Anderson, Eric Polley, Larry Rubinstein, Apurva Srivastava, Deborah Wilsker, et al.  The National Cancer Institute ALMANAC: a comprehensive screening resource for the detection of anticancer drug pairs with enhanced therapeutic activity.  Cancer Research, 77(13):3564–3576, 2017.
  * [25] Jennifer O’Neil, Yair Benita, Igor Feldman, Melissa Chenard, Brian Roberts, Yaping Liu, Jing Li, Astrid Kral, Serguei Lejnine, Andrey Loboda, et al.  An unbiased oncology compound screen to identify novel combination strategies.  Molecular Cancer Therapeutics, 15(6):1155–1162, 2016.
  * [26] Harry L Morgan.  The generation of a unique machine description for chemical structures-a technique developed at chemical abstracts service.  Journal of Chemical Documentation, 5(2):107–113, 1965.
  * [27] Ethan Perez, Florian Strub, Harm De Vries, Vincent Dumoulin, and Aaron Courville.  Film: Visual reasoning with a general conditioning layer.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 32, 2018.
  * [28] Hiral A Shah, James H Fischer, Neeta K Venepalli, Oana C Danciu, Sonia Christian, Meredith J Russell, Li C Liu, James P Zacny, and Arkadiusz Z Dudek.  Phase i study of aurora a kinase inhibitor alisertib (mln8237) in combination with selective vegfr inhibitor pazopanib for therapy of advanced solid tumors.  American journal of clinical oncology, 42(5):413–420, 2019.
  * [29] Shaoming Zhu, Tian Zhang, Lei Zheng, Hongtao Liu, Wenru Song, Delong Liu, Zihai Li, and Chong-xian Pan.  Combination strategies to maximize the benefits of cancer immunotherapy.  Journal of Hematology & Oncology, 14(1):1–33, 2021.
  * [30] Qi Feng, Aoshuang Zhou, Huachun Zou, Suzanne Ingle, Margaret T May, Weiping Cai, Chien-Yu Cheng, Zuyao Yang, and Jinling Tang.  Quadruple versus triple combination antiretroviral therapies for treatment naive people with hiv: systematic review and meta-analysis of randomised controlled trials.  bmj, 366, 2019.
  * [31] Cordula Hirsch and Stefan Schildknecht.  In vitro research reproducibility: Keeping up high standards.  Frontiers in pharmacology, 10:1484, 2019.
  * [32] Ziga Avsec, Vikram Agarwal, Daniel Visentin, Joseph R Ledsam, Agnieszka Grabska-Barwinska, Kyle R Taylor, Yannis Assael, John Jumper, Pushmeet Kohli, and David R Kelley.  Effective gene expression prediction from sequence by integrating long-range interactions.  bioRxiv, 2021.
  * [33] Xin Shao, Ning Lv, Jie Liao, Jinbo Long, Rui Xue, Ni Ai, Donghang Xu, and Xiaohui Fan.  Copy number variation is highly correlated with differential gene expression: a pan-cancer study.  BMC medical genetics, 20(1):1–14, 2019.
  * [34] Jake P Taylor-King.  Rethinking rare disease: longevity-enhancing drug targets through x-linked aneuploidy.  Trends in Genetics, 2021.
  * [35] Thomas Gaudelet, Ben Day, Arian R Jamasb, Jyothish Soman, Cristian Regep, Gertrude Liu, Jeremy BR Hayter, Richard Vickers, Charles Roberts, Jian Tang, et al.  Utilizing graph machine learning within drug discovery and development.  Briefings in bioinformatics, 22(6):bbab159, 2021.
  * [36] Aravind Subramanian, Rajiv Narayan, Steven M Corsello, David D Peck, Ted E Natoli, Xiaodong Lu, Joshua Gould, John F Davis, Andrew A Tubelli, Jacob K Asiedu, et al.  A next generation connectivity map: L1000 platform and the first 1,000,000 profiles.  Cell, 171(6):1437–1452, 2017.
  * [37] Haide Chen, Yuan Liao, Guodong Zhang, Zhongyi Sun, Lei Yang, Xing Fang, Huiyu Sun, Lifeng Ma, Yuting Fu, Jingyu Li, et al.  High-throughput microwell-seq 2.0 profiles massively multiplexed chemical perturbation.  Cell discovery, 7(1):1–4, 2021.
  * [38] Stefan Peidli, Tessa Durakis Green, Ciyue Shen, Torsten Gross, Joseph Min, Jake Taylor-King, Debora Marks, Augustin Luna, Nils Bluthgen, and Chris Sander.  scperturb: Information resource for harmonized single-cell perturbation data.  bioRxiv, 2022.
  * [39] Mark-Anthony Bray, Shantanu Singh, Han Han, Chadwick T Davis, Blake Borgeson, Cathy Hartland, Maria Kost-Alimova, Sigrun M Gustafsdottir, Christopher C Gibson, and Anne E Carpenter.  Cell painting, a high-content image-based assay for morphological profiling using multiplexed fluorescent dyes.  Nature protocols, 11(9):1757–1774, 2016.
  * [40] Tuan-Minh Nguyen, Adib Shafi, Tin Nguyen, and Sorin Draghici.  Identifying significantly impacted pathways: a comprehensive review and assessment.  Genome biology, 20(1):1–15, 2019.
  * [41] Jake P Taylor-King, Etienne Baratchart, Andrew Dhawan, Elizabeth A Coker, Inga Hansine Rye, Hege Russnes, S Jon Chapman, David Basanta, and Andriy Marusyk.  Simulated ablation for detection of cells impacting paracrine signalling in histology analysis.  Mathematical medicine and biology: a journal of the IMA, 36(1):93–112, 2019.
  * [42] Patricia Jaaks, Elizabeth A Coker, Daniel J Vis, Olivia Edwards, Emma F Carpenter, Simonetta M Leto, Lisa Dwane, Francesco Sassi, Howard Lightfoot, Syd Barthorpe, et al.  Effective drug combinations in breast, colon and pancreatic cancer cells.  Nature, 603(7899):166–173, 2022.
  * [43] Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell.  Simple and scalable predictive uncertainty estimation using deep ensembles, 2017.
  * [44] Moksh Jain, Salem Lahlou, Hadi Nekoei, Victor Butoi, Paul Bertin, Jarrid Rector-Brooks, Maksym Korablyov, and Yoshua Bengio.  DEUP: direct epistemic uncertainty prediction.  CoRR, abs/2102.08501, 2021.
  * [45] Olivier Borkowski, Mathilde Koch, Agnès Zettor, Amir Pandi, Angelo Cardoso Batista, Paul Soudier, and Jean-Loup Faulon.  Large scale active-learning-guided exploration for in vitro protein production optimization.  Nature Communications, 11(1):1872, Dec 2020.
  * [46] Ross D. King, Kenneth E. Whelan, Ffion M. Jones, Philip G. K. Reiser, Christopher H. Bryant, Stephen H. Muggleton, Douglas B. Kell, and Stephen G. Oliver.  Functional genomic hypothesis generation and experimentation by a robot scientist.  Nature, 427(6971):247–252, Jan 2004.
  * [47] Pablo Carbonell, Adrian J. Jervis, Christopher J. Robinson, Cunyu Yan, Mark Dunstan, Neil Swainston, Maria Vinaixa, Katherine A. Hollywood, Andrew Currin, Nicholas J. W. Rattray, Sandra Taylor, Reynard Spiess, Rehana Sung, Alan R. Williams, Donal Fellows, Natalie J. Stanford, Paul Mulherin, Rosalind Le Feuvre, Perdita Barran, Royston Goodacre, Nicholas J. Turner, Carole Goble, George Guoqiang Chen, Douglas B. Kell, Jason Micklefield, Rainer Breitling, Eriko Takano, Jean-Loup Faulon, and Nigel S. Scrutton.  An automated design-build-test-learn pipeline for enhanced microbial production of fine chemicals.  Communications Biology, 1(1):66, Dec 2018.
  * [48] Brian Hie, Bryan D Bryson, and Bonnie Berger.  Leveraging uncertainty in machine learning accelerates biological discovery and design.  Cell Systems, 11(5):461–477, 2020.
  * [49] Mark Davies, Michał Nowotka, George Papadatos, Nathan Dedman, Anna Gaulton, Francis Atkinson, Louisa Bellis, and John P Overington.  ChEMBL web services: streamlining access to drug discovery data and utilities.  Nucleic Acids Research, 43(W1):W612–W620, 2015.
  * [50] Sunghwan Kim, Jie Chen, Tiejun Cheng, Asta Gindulyte, Jia He, Siqian He, Qingliang Li, Benjamin A Shoemaker, Paul A Thiessen, Bo Yu, et al.  Pubchem 2019 update: improved access to chemical data.  Nucleic acids research, 47(D1):D1102–D1109, 2019.
  * [51] Mahmoud Ghandi, Franklin W. Huang, Judit Jané-Valbuena, Gregory V. Kryukov, Christopher C. Lo, E. Robert McDonald, Jordi Barretina, Ellen T. Gelfand, Craig M. Bielski, Haoxin Li, Kevin Hu, Alexander Y. Andreev-Drakhlin, Jaegil Kim, Julian M. Hess, Brian J. Haas, François Aguet, Barbara A. Weir, Michael V. Rothberg, Brenton R. Paolella, Michael S. Lawrence, Rehan Akbani, Yiling Lu, Hong L. Tiv, Prafulla C. Gokhale, Antoine de Weck, Ali Amin Mansour, Coyin Oh, Juliann Shih, Kevin Hadi, Yanay Rosen, Jonathan Bistline, Kavitha Venkatesan, Anupama Reddy, Dmitriy Sonkin, Manway Liu, Joseph Lehar, Joshua M. Korn, Dale A. Porter, Michael D. Jones, Javad Golji, Giordano Caponigro, Jordan E. Taylor, Caitlin M. Dunning, Amanda L. Creech, Allison C. Warren, James M. McFarland, Mahdi Zamanighomi, Audrey Kauffmann, Nicolas Stransky, Marcin Imielinski, Yosef E. Maruvka, Andrew D. Cherniack, Aviad Tsherniak, Francisca Vazquez, Jacob D. Jaffe, Andrew A. Lane, David M. Weinstock, Cory M. Johannessen, Michael P. Morrissey, Frank Stegmeier, Robert Schlegel, William C. Hahn, Gad Getz, Gordon B. Mills, Jesse S. Boehm, Todd R. Golub, Levi A. Garraway, and William R. Sellers.  Next-generation characterization of the Cancer Cell Line Encyclopedia.  Nature, 569(7757):503–508, May 2019.
  * [52] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.
  * [53] RDKit: Open-source cheminformatics.  http://www.rdkit.org.

##  Appendix A Methods

###  A.1. Model description

We frame the problem of pairwise drug synergy prediction as a regression task
$(\\{a,b\\},\hat{\textnormal{s}})$: given a pair of drugs $a,b$, we aim to
predict their (pooled) level of synergy, $\hat{\textnormal{s}}$. Our proposed
architecture is an end-to-end framework trained with a mean square error (MSE)
criterion.

Our model can be decomposed into two modules. First, a single drug module,
$E$, which produces representations (or embeddings) for the drugs based on
their chemical structure information. The embeddings from a pair of drugs are
used as input to the combination module $P$, which directly estimates the
synergy score; see Figure S1.

Figure S1. Overview of the RECOVER model. Drug representations are fed into
the Single drug module which is composed of an MLP which can be conditioned on
cell line. Given two drug embeddings, the synergy is predicted using the
Combination module, composed of a bilinear operation followed by an MLP. The
Combination MLP can be conditioned on cell line as well.

Further, uncertainty estimation methods are used in order to estimate the
predictive distribution of synergies $p(\hat{\textnormal{s}}|\\{a,b\\})$ for
each drug pair $\\{a,b\\}$, as opposed to a point estimate. The predictive
distributions of drug pairs are given as input to an acquisition function in
order to decide which combinations should be tested in vitro, balancing
between combinations that are informative, i.e. that can reduce the
generalization error of the model later on, and combinations that are likely
to be synergistic.

####  A.1.1. Single drug module

Let $X_{\mathcal{D}}\in\mathbb{R}^{n_{\mathcal{D}}\times l_{\mathcal{D}}}$
denote the matrix of drug features, where $n_{\mathcal{D}}$ is the number of
drugs in $\mathcal{D}$ and $l_{\mathcal{D}}$ corresponds to the number of raw
features that describe each drug. Drug features used in this work include
molecular fingerprints [10] and one-hot encoding of the drugs.

The single drug module can be written as a function
$E\mathrel{\mathop{\mathchar
58\relax}}\mathcal{D}\rightarrow\mathbb{R}^{k_{\mathcal{D}}}$ where
$k_{\mathcal{D}}$ corresponds to the dimension of the output vector
representation (or embedding) of each drug. Our single drug module is a simple
multi-layer perceptron (MLP) that takes raw features of drugs as input and
outputs an updated vector representation. This MLP can be conditioned on cell
line as described in Appendix A.1.3.

####  A.1.2. Combination module

Given a set of drugs $\mathcal{D}$, the combination module corresponds to a
function $P\mathrel{\mathop{\mathchar
58\relax}}\binom{\mathcal{D}}{2}\mapsto\mathbb{R}$ that maps a pair of drugs
to their Bliss synergy score. We remark first that $P$ should be agnostic to
the order of the two drugs. Hence, the first operation of $P$ correspond to a
permutation invariant function – such as element-wise sum, mean, or max
operations – applied to the two vector representations corresponding to each
drug. In this work, we use a bilinear operation defined by a tensor
$B\in\mathbb{R}^{k_{\mathcal{D}}\times k_{\mathcal{D}}\times k}$, where $k$ is
a hyperparameter corresponding to the dimension of the vector representation
of a drug combination. To ensure permutation invariance, we enforce that every
slice across the third dimension (denoted as $B_{i}$), is a symmetric matrix.
Note that we do not enforce $B_{i}$ to be positive definite, hence $B_{i}$
does not necessarily define a scalar product. The output of this permutation
invariant function is fed to an MLP that outputs the predicted synergy for the
pair of drugs. Again, the MLP can be conditioned on cell lines, see Appendix
A.1.3.

####  A.1.3. Cell line conditioning

As a drug effect is context dependent, the synergy of a combination of two
drugs can be different in experiments using different cell lines. To account
for the cell line in our model we condition upon it using FiLM [27]. In
essence, the FiLM approach learns an affine transformation of the activation
of each neuron in the MLP.

We denote the matrix of cell line features by
$X_{\mathcal{C}}\in\mathbb{R}^{n_{\mathcal{C}}\times l_{\mathcal{C}}}$, with
$\mathcal{C}$ the set of cell lines, $n_{\mathcal{C}}$ corresponding to the
number of cell lines in $\mathcal{C}$ and $l_{\mathcal{C}}$ giving the number
of raw features for each cell line.

The feature representation of the cell line is either based on a one-hot
encoding, or on information about mutations and basal level of gene
expression. The former approach relies on having data for each cell line in
the training set and cannot generalise to new cell lines. The second approach
makes use of features that represent cell lines as described in Appendix A.4.

###  A.2. Searching the space of drug combinations

####  A.2.1. Uncertainty estimation

Estimating the uncertainty of the predictions is a key step towards providing
reliable recommendations as well as driving the exploration with SMO. For this
purpose, we use a common uncertainty estimation method: deep ensembles [43].
Given an ensemble of models which differ only in the initialization of the
parameters, the predictions of the different models are considered as samples
from the predictive distribution. In this work, we define uncertainty as the
standard deviation of the predictive distribution, and can be estimated from
the standard deviation between the predictions of the different members of the
ensemble. Unless specified otherwise, we use a deep ensemble of size 5 as the
uncertainty estimation method in our in silico experiments, and of size 36 for
recommendation generation, as described in Appendix A.3.

Note that for completeness, we investigated other methods for uncertainty
quantification in some of the in silico experiments, including direct
estimation of the standard deviation of the predictive distribution — in a
similar fashion to Direct Epistemic Uncertainty Prediction (DEUP) [44], see
Appendix C.1 for more details.

####  A.2.2. Sequential model optimization

Sequential model optimization (SMO) aims at discovering an input
$x^{\star}\in\mathcal{X}$ maximizing an objective function $S$:

(3) $x^{\star}\in\operatorname*{arg\,max}\limits_{x\in\mathcal{X}}S(x).$

The SMO approach consists in tackling this problem by iteratively querying the
objective function $S$ in order to find a maximizer $x^{\star}$ in a minimal
number of steps. At each step $t$, the dataset is augmented such that
$\mathcal{D}_{t}$ contains all the inputs that have already been acquired at
time $t$. The dataset $\mathcal{D}_{t}$ is then used to find the next query
$x^{(t+1)}$. In the context of drug combinations, $x$ corresponds to a pair of
drugs, and the objective function $S$ corresponds to the synergy score.

SMO has been prospectively applied to: optimize the production of proteins in
cell free systems [45]; determine gene functions in yeast [46]; enhance the
production of fine chemicals in Escherichia coli [47]; and to identify
inhibitors of Mycobacterium tuberculosis growth [48].

In what follows, $f$ refers to an estimator of the objective function $S$. One
may notice that several properties of the potential queries $x^{(t)}$ should
be taken into account. One would like to find an $x^{(t)}$ that would be
informative to acquire (i.e., the uncertainty at $x^{(t)}$ is high) in order
to obtain a reliable estimator of the objective function early on. On the
other hand, one would like to find an $x^{(t)}$ that is a good guess in the
sense that $f(x^{(t)})$ is close to the expected maximum
$\max_{x\in\mathcal{X}}f(x)$. Looking for queries which are informative is
referred to as exploration while looking for queries which are expected to
maximize the objective function is called exploitation.

The key challenge of SMO is to balance between exploration and exploitation.
This is typically achieved by designing an acquisition function (or strategy)
$\alpha$ which defines a score on the space of inputs $\mathcal{X}$ and takes
into account both the expected $f(x)$ and an estimate of the uncertainty at
$x$. The input which maximizes the score $\alpha$ is chosen as the next query.
An overview of the SMO approach is presented in Algorithm 1.

Input: Initial data $\mathcal{D}_{0}$, objective function estimator $f$

for $t\in\\{1,2,...\\}$ do

Select new $x^{(t+1)}$ by optimizing an acquisition function $\alpha$
$x^{(t+1)}=\operatorname*{arg\,max}\limits_{x}\alpha(x;f)$

Query objective function $S$ to obtain $y^{t+1}$

Augment data $\mathcal{D}_{n+1}=\mathcal{D}_{n}\cup\\{(x^{t+1},y^{t+1})\\}$

Update estimator $f$

end for

Algorithm 1 Sequential model optimization

In what follows, we assume that we have access to an estimate of the mean of
the predictive distribution, $\hat{\mu}(x)$, as well as an estimate of the
uncertainty $\hat{\sigma}(x)$. The key acquisition functions considered are
detailed below.

##### Brute-force

$\alpha(x)$ corresponds to random noise, and therefore the drug combinations
are selected at random.

##### Greedy

$\alpha(x)=\hat{\mu}(x)$. This acquisition function corresponds to pure
exploitation whereby we select drug combinations with the highest predicted
synergy.

##### Pure exploration

$\alpha(x)=\hat{\sigma}(x)$. This acquisition function corresponds to pure
exploration. The strategy aims at labelling the most informative examples in
order to reduce model uncertainty as fast as possible, and corresponds to the
traditional strategy in Active Learning.

##### Upper confidence bound (UCB)

$\alpha(x)=\mu(x)+\kappa\hat{\sigma}(x)$. This strategy balances between
exploration and exploitation. $\kappa\in\mathbb{R}$ is a hyperparameter that
is typically positive. Higher values of $\kappa$ give more importance to
exploration.

Unless specified otherwise, in silico experiments involving SMO were performed
using UCB with $\kappa=1$.

In our synthetic SMO experiments, the model was reinitialized and trained from
scratch on all visible data after each query. Whilst not designed for optimal
computational efficiency, this procedure ensures that the model is not
overfitting on examples that have been acquired early on.

###  A.3. Recommendation generation

In order to generate the recommendations for in vitro experiments, we trained
3 models using 3 different seeds on the NCI-ALMANAC study, restricting
ourselves to samples from the MCF7 cell line. We refer to these 3 models as
pretrained.

Afterwards, we fine-tune using prospectively generated data. More precisely,
the weights of one of the pretrained models were loaded, and some additional
training was performed on prospectively generated data only, using early
stopping. This fine-tuning process was repeated with 12 different seeds for
each pretrained model. The end result being that we obtain an ensemble of 36
fine-tuned models in total.

This ensemble was used to generate predictions $(\hat{\mu},\hat{\sigma})$ for
all candidate combinations. We then use UCB with $\kappa=1$ to obtain a score
according to which all candidates were ranked.

###  A.4. Dataset processing

Datasets have been pre-processed and normalized in a centralised data
repository RESERVOIR111see http:

github.com
RECOVERcoalition
RESERVOIR. The repository unifies data around biologically relevant molecules
and their interactions. Pre-processing and normalizing scripts are provided
for traceability, and a Python API has been made available to facilitate
access. Below, the major data types included are briefly described.

####  A.4.1. Drugs

Data on drugs and biologically active compounds has been extracted from Chembl
[49], pre-processed and indexed with unique identifiers. A translation engine
has been provided such that a compound can be translated to a unique
identifier using generic or brand drug names, SMILES strings and Pubchem [50]
CIDs.

####  A.4.2. Cell Line Features

Additionally, RESERVOIR retrieved cell line features from the Cancer
Dependency Map [51]. These include genetic mutations, base level gene
expression and metadata.

####  A.4.3. Drug combinations

Literature drug combination data was extracted from DrugComb version 1.5 [22].
Quality control was applied on the experiments in DrugComb. Only blocks (i.e.
combination matrices) complying to the following criteria were selected: (a.)
filter out erroneous blocks that show very low variance, specifically
inhibition standard deviation $\leq$0.05, (b.) filter out small blocks less
than 3 $\times$ 3 dimensions, (c.) filter out blocks with extreme inhibition
values, such that 5% $<$ mean pooled growth inhibition $<$95%.

The dataset used for model pretraining and in silico experiments consists of
4463 data points relative to experiments on MCF7 cell line expressed as max
Bliss which were reported in the Almanac study. These data correspond to 4271
unique drug combinations made up by 95 unique drugs.

The prediction set for experiment selection was built by taking 54 out of the
95 Almanac drugs for which a mode of action (MoA) was annotated in ChEMBL 25
[49]. An additional 54 drugs were obtained by clustering 719 drugs with known
MoA that are included in DrugComb but are not part of Almanac. Clustering was
performed with the $k$-medoids algorithm as implemented in scikit-learn 0.24.2
[52] (n_clusters=54, metric=Tanimoto similarity, init=k-medoids++), drugs were
encoded by Morgan fingerprints with radius 2 and 1024 bits calculated with
RDKit [53]. A representative compound for each cluster was obtained by taking
the cluster centroid.

Three of the centroid drugs were replaced due to lack of availability from
commercial vendors or due to poor reported solubility. Replacements for each
of the three drugs were selected by taking the nearest analogue (evaluated by
Tanimoto similarity) in the same cluster. 54 Almanac and 54 non-Almanac
compounds thus selected were used to build a set of 2916 binary combinations
made up by one Almanac and one non-Almanac compound.

###  A.5. Experimental protocol

####  A.5.1. Cell lines

MCF-7 cells were obtained from ATCC and maintained in DMEM (Corning)
supplemented with 10% FBS (Corning) and Antibiotic-Antimycotic (Gibco) at 37°C
in 5% $\text{CO}_{2}$ in a humidified incubator. Before the screens, the cell
lines were passaged twice after thawing. Cultures were confirmed to be free of
mycoplasma infection using the MycoAlert Mycoplasma Detection Kit (Lonza).

####  A.5.2. Drug Combination Assays

A list of the tested compounds along with their concentration ranges and
target mechanisms is given as a supplementary table. The compound-specific
concentration ranges were selected based on their published activities. All
compounds were pre-diluted in DMSO to a stock concentration that varied from
10 to 50 mM, depending on the final concentration range required for each
compound.

Briefly, compounds were plated as a 6 $\times$ 6 dose-response combination
matrix in natural 384 well plates (Greiner), in a serial 1:3 dilutions of each
agent (5 concentrations) and only DMSO as the lowest concentration. We used a
combination plate layout where six compound pairs could be accommodated on one
384 well plate. A set of control wells with DMSO was included on all plates as
negative control. To ensure reproducibility and comparability with the
subsequent combination studies, the IC50 of Doxorubicin was used as reference
in a 6-point dose response format in each plate as positive (total killing)
control.

Cells were seeded in white 384-well plates (Greiner) at 1000 cells/well in 50
$\mu$L of media using a multidrop dispenser and allowed to attach for 2 h.
Compounds from pre-plated matrix plates were transferred to each well using a
100 nL head affixed to an Agilent Bravo automated liquid handling platform,
and plates were incubated at 37°C in 5% CO2 for an additional 72 hours. To
measure the cell viability, CellTiter-Glo reagent (diluted 1:6 in water,
Promega) was dispensed into the wells (30 $\mu$L), incubated for 3 minutes,
and luminescence was read on a Envision plate reader (Perkin-Elmer). Final
DMSO concentration in assay wells was 0.2%. The assay was performed with 3
biological replicates.

Supplementary Figures

\begin{overpic}[width=433.62pt]{img/appendix/FigureS2_ALMANAC_dimensions.pdf}
\put(5.0,65.0){A}
\end{overpic}\begin{overpic}[width=433.62pt]{img/appendix/FigureS2_ALMANAC_drug_phases.pdf}
\put(5.0,80.0){B} \end{overpic}

\begin{overpic}[width=433.62pt]{img/appendix/FigureS2_ALMANAC_synergy_matrix.pdf}
\put(0.0,82.0){C} \end{overpic}

\begin{overpic}[width=433.62pt]{img/appendix/FigureS2_ALMANAC_drug_synergies.pdf}
\put(-5.0,93.5){D}
\end{overpic}\begin{overpic}[width=433.62pt]{img/appendix/FigureS2_ONEIL2ALMANAC_corr.pdf}
\put(-5.0,90.0){E} \end{overpic}

Figure S2. Model performance was evaluated on synergy data from the NCI-
ALMANAC viability screen [24] for the MCF7 cell line. After quality control,
synergy scores could be computed for 4271 unique drug pairs. (A.) Distribution
of dose-response matrix dimensions. (B.) Distribution of clinical phases of
the considered 95 distinct drugs. (C.) Synergy scores for each drug-pair.
White squares indicate missing or removed (low-quality) dose-response
matrices. (D.) Histogram of average synergy scores for each drug computed
across all drug pairs recorded. (E.) Scatter plot across 783 (drug-pair, cell
line) tuples found in both NCI-ALMANAC and O’Neil 2016 [25].
\begin{overpic}[width=390.25534pt]{img/appendix/FigureS3.png} \end{overpic}
Figure S3. Bliss synergy dose response plots the top 12 synergistic drug
combinations (excluding drug combos found in Figure S4).
\begin{overpic}[width=390.25534pt]{img/appendix/FigureS4.png}
\put(2.0,38.0){A} \put(50.5,38.0){B} \end{overpic} Figure S4. Bliss synergy
dose response plots for two drug combinations discovered within the SMO Search
1 round: (A.) Alisertib & Pazopanib and (B.) Flumatinib & Mitoxantrone.

Supplementary Information

##  Appendix B Model development and evaluation

We investigated various aspects of the performance of RECOVER for the
prediction of Bliss synergy scores. All results presented in this section have
been computed on the NCI-ALMANAC study restricted to the MCF7 cell line.
Combinations are split randomly into training/validation/test (70%/20%/10%).
We restrict ourselves to MCF7 for consistency with prospective in vitro
experiments.

###  B.1. RECOVER outperforms baseline models

We first compare our approach to several baseline models. In Table 2, we
demonstrate that RECOVER outperforms baseline models in terms of $R^{2}$ and
Spearman correlation metrics. For this evaluation, the hyperparameters of
RECOVER have been optimized within the following set of values. Because it was
not tractable to perform a grid search over all possible values at once,
hyperparameters have been optimized one at a time in an iterative way. The set
of parameters that yielded best performance is highlighted, and used for all
following experiments (both in silico and in vitro experiments):

  * •

Learning rate: [$1\times 10^{-1}$, $1\times 10^{-2}$, $1\times 10^{-3}$,
$\mathbf{1\times 10^{-4}}$, $1\times 10^{-5}$, $1\times 10^{-6}$]

  * •

Batch size: [16, 32, 64, 128, 256]

  * •

Weight decay: [$1$, $1\times 10^{-1}$, $\mathbf{1\times 10^{-2}}$, $1\times
10^{-3}$, $1\times 10^{-4}$, $1\times 10^{-5}$, $1\times 10^{-6}$, $0$]

  * •

Morgan fingerprint radius: [2, 3, 4, 5]

  * •

Morgan fingerprint dimension: [1024, 2048]

  * •

Output dimension of the single drug MLP: [16, 32, 64, 128, 256]

  * •

Dimension(s) of the hidden layer(s) of the single drug MLP: [[512], [256],
[1024], [2048], [4096], [1024, 1024], [1024, 1024, 1024], [1024, 512], [1024,
512, 256]]

  * •

Dimension(s) of the hidden layer(s) of the combination MLP: [[32], [64],
[128], [256], [64, 16], [64, 32], [64, 64]]

Moreover, a grid search was performed to optimize the hyperparameters of the
Gradient Boosting Trees baseline model. The number of trees was set to $100$.
The set of parameters that yielded best performance is highlighted:

  * •

Maximum tree depth: [2, 5, 10, 20]

  * •

Minimum number of samples to split a node: [2, 5, 10, 20, 50]

  * •

Learning rate: [0.0001, 0.001, 0.01, 0.1, 1]

  * •

Maximum feature: [all, $\sqrt{\text{total number of features}}$,
$\log_{2}[\text{total number of features}]$]

Similarly, a grid search was performed for the Linear SVM baseline model:

  * •

Tolerance for stopping criterion: [$1\times 10^{-1}$, $1\times 10^{-2}$,
$\mathbf{1\times 10^{-3}}$, $1\times 10^{-4}$, $1\times 10^{-5}$]

  * •

Regularization C: [0.0001, 0.001, 0.01, 0.1, 1., 10]

| $R^{2}$ | Spearman corr.
---|---|---
RECOVER | $\mathbf{0.242\pm 0.006}$ | $\mathbf{0.466\pm 0.007}$
Linear SVM | $0.171\pm 0.01$ | $0.453\pm 0.004$
Gradient Boosting Trees | $0.172\pm 0.002$ | $0.459\pm 0.001$
Table 2. Comparison of RECOVER with baseline models for the synergy score
prediction within the MCF7 cell line. Standard deviation computed over 3
seeds.

###  B.2. Feature importance study

Through investigation of different drug features, we find a large proportion
of the performance of RECOVER can be achieved given the identity of the drugs
alone, and that structural information allows for a slight increase in
performance. As shown in Table 3, the performance of the model is similar
whether the one hot encoding of the drug or its Morgan fingerprint is used as
input. We notice a slight improvement when using both feature types together.
Note that the number of parameters of the model is always the same regardless
of the type of feature provided as input. When a feature type is not used, the
corresponding part of the drug feature vector is set to zero without changing
the underlying dimension.

| $R^{2}$ | Spearman corr.
---|---|---
RECOVER (fingerprint + one hot) | $\mathbf{0.242\pm 0.006}$ | $\mathbf{0.466\pm 0.007}$
RECOVER (fingerprint) | $0.232\pm 0.007$ | $0.458\pm 0.010$
RECOVER (one hot) | $0.230\pm 0.029$ | $0.449\pm 0.017$
Table 3. Feature importance study for the prediction of max Bliss score on the
MCF7 cell line. Standard deviation computed over 3 seeds.

###  B.3. Upper bounds on model performance

We investigate RECOVER performance with regards to Spearman correlation and
$R^{2}$. Whilst predictive power appears modest, we are still able to identify
highly synergistic drug combinations in simulated SMO experiments, see
Appendix B.4. Several aspects that may limit predictive power: experimental
noise, and nonuniformity of maximum Bliss synergy scores.

In Figure S2C, we note most data points are close to zero, with some examples
very far from the mean, i.e., the examples of interest. As an example, let us
consider the case of Spearman correlation. Given that the observations are
noisy, the observed rank among synergies might be corrupted compared to the
true ordering — especially in the region close to zero where the density of
examples is very high.

The non-uniformity of synergy scores leads to some difficulties in evaluating
fairly the performance of RECOVER. For example, the positive tail of the
distribution, which is the region of interest, represent a very small
percentage of the total number of examples and thus have a little effect on
the value of the aggregated statistic.

In order to get a better understanding of the performance of our model, we
compare the reported aggregated statistics to an upper bound which takes into
account the presence of noise in the observations in addition to the
distribution of synergy scores.

##### Evaluating the level of noise

Consider all replicates from the NCI-ALMANAC study. Two examples are
considered replicates when the same pair of drugs has been tested on the same
cell line. We found $1960$ triplets ({drug1, drug2}, cell line) that had been
tested several times. For each triplet, we computed the standard deviation of
the maximum Bliss score across the replicates. We refer to this as the level
of noise for a given triplet. We then computed the average level of noise
$\bar{\eta}$ across all triplets.

##### Estimating upper bounds

Given an average level of noise $\eta$ and the distribution of synergy scores
in NCI-ALMANAC, we simulated a noisy acquisition process as follows: the
synergies from NCI-ALMANAC were considered as the true synergies, and noisy
observations were obtained by corrupting the true synergies with some Gaussian
noise $\mathcal{N}(0,\eta^{2})$. We then considered a perfect regression model
which fits the noisy observations exactly, and evaluated its performance on
the true synergies. Upper bounds are defined by the performance of this
perfect regression model.

Upper bounds have been computed for $R^{2}$ and Spearman correlation using
various levels of noise and are reported in Figure S5. We see that the noisy
acquisition process alone leads to significant limitations in the performance
that can be reached. While there is still room for improvement, the
performance of RECOVER is reasonably close to the hypothetical maximum. For
example, RECOVER achieves 0.47 Spearman correlation, while the highest
achievable Spearman correlation is estimated to be 0.64.

\begin{overpic}[width=195.12767pt]{img/appendix/FigureS5_up_boun_r2.png}
\put(2.0,65.0){A} \end{overpic}

\begin{overpic}[width=195.12767pt]{img/appendix/FigureS5_up_boun_spear.png}
\put(2.0,65.0){B} \end{overpic}

Figure S5. Comparison of the performance of RECOVER with upper bounds that
take into account the level of noise in the data, and the non uniform
distribution of synergies: (A.) $R^{2}$ value, (B.) Spearman correlation. The
upper bound, shown in blue, is a function of the level of noise. The standard
deviation for the upper bound is lower than $10^{-3}$ for both statistics
(estimated over 3 seeds). The estimated level of noise in our dataset is shown
in green. The performance of RECOVER as well as its standard error are shown
in red.

###  B.4. Performance analysis on the tail of the distribution of synergies

We now show that RECOVER achieves good performance on the positive tail of the
distribution of synergies, which is necessary to successfully identify highly
synergistic combinations within a SMO setting.

Given a model trained on NCI-ALMANAC, we query the top $k$ combinations with
highest predicted synergy within the test set, and compute the percentage of
combinations which are truly synergistic (arbitrarily defined by a maximum
Bliss synergy score of above $30$) within queried examples. We refer to $k$ as
the size of the query. In Figure S6, we report the percentage of synergistic
combinations as a function of the size of the query. The percentage of
synergistic combinations is superior to the proportion in the whole test set,
meaning that the model performs far beyond the level of randomness. For
instance, with a query size of $30$, we observe a $\sim$5-fold enrichment in
synergistic combinations. However, this enrichment decreases with the size of
the query.

Figure S6. Percentage of highly synergistic combinations (synergy score $>30$)
in queried combinations using RECOVER as a function of the size of the query.
Shaded area corresponds to the standard deviation computed over 3 seeds.

##  Appendix C In silico Sequential Model optimization

We benchmarked the SMO pipelines, whereby the model is shown a fraction of the
full data set and can choose sample points to unblind. Our in silico
experiments try to mirror as closely as possible the setting of the in vitro
experiments. Therefore, unless specified otherwise, experiments are restricted
to the MCF7 cell line and 30 combinations were acquired at a time using the
same model as the one used to generate recommendations for the in vitro
experiments. Uncertainty is estimated using a deep ensemble of size 5, unless
stated otherwise.

###  C.1. Alternative uncertainty estimators: Direct Estimation

As deep ensembles are only one method for uncertainty quantification, we
tested another approach: directly estimating the uncertainty in the style of
DEUP [44]. Here, two models are initialized: the first one, denoted as the
mean predictor, predicts the expected synergy $\hat{\mu}$ and is trained with
Mean Square Error (MSE). The second model, denoted as the uncertainty
predictor, outputs an estimate of the standard deviation of the predictive
distribution $\hat{\sigma}$ and is trained to minimize the following negative
log-likelihood criterion:

(4)
$NLL=\frac{\log(\hat{\sigma}^{2})}{2}+\frac{(y-\hat{\mu})^{2}}{2\hat{\sigma}^{2}},$

where $y$ is the ground truth variable of interest. This criterion allows us
to get an estimator $\hat{\sigma}$ of the standard deviation of the predictive
distribution for each combination: given the expected synergy $\hat{\mu}$
predicted by the mean predictor, and the actual observation $y$, we wish to
find $\hat{\sigma}$ that maximizes the probability of $y$ assuming a
predictive distribution of the form $\mathcal{N}(\hat{\mu},\hat{\sigma})$.

(5) $\displaystyle\log p(y|\hat{\mu},\hat{\sigma})$
$\displaystyle=\log\left[\frac{1}{\sqrt{2\pi\hat{\sigma}^{2}}}e^{-\frac{1}{2}(\frac{y-\hat{\mu}}{\hat{\sigma}})^{2}}\right]$
(6)
$\displaystyle=-\frac{1}{2}\log(2\pi)-\frac{1}{2}\log(\hat{\sigma}^{2})-\frac{1}{2}\frac{(y-\hat{\mu})^{2}}{\hat{\sigma}^{2}}$

Therefore, maximizing $\log p(y|\hat{\mu},\hat{\sigma})$ w.r.t. $\hat{\sigma}$
is equivalent to minimizing the $NLL$ criterion presented in Equation 4. Note
that only the uncertainty predictor is trained using this criterion. The mean
predictor is trained using Mean Square Error (MSE) as we found experimentally
that it was more stable. When fixing $\hat{\sigma}=1$, $NLL$ corresponds to
the MSE criterion.

###  C.2. Backtesting RECOVER demonstrates efficient exploration of drug
combination space

In order to simulate real-world interactions with the wet lab, we start with a
set of 30 randomly chosen drug pairs that the pipeline is initially trained
on, while the rest of the data is hidden from the pipeline. Hence, these
combinations are always part of the visible set. We now perform an iteration:
we split the visible set into a training and validation set (80/20), train the
model using early stopping and then acquire 30 additional drug pairs from the
hidden set. The entire procedure is repeated again with the new training set
of $60$ drug pairs (i.e., $30+30$), and then with $90$, and so on, until no
more drug pairs are left to acquire from the hidden set. A range of methods
can be used to generate recommendations for new drug pairs, typically as
specified by an acquisition function [23]. The acquisition functions that we
investigate are detailed in Appendix A.2.2.

We compare acquisition functions by the rate at which they unblind the set of
top 1% of synergistic combinations in the NCI-ALMANAC dataset within each
synthetic experimental round (or iteration). As shown in Figure S7A, Greedy
acquisition and Upper Confidence Bound (UCB) acquisition, two model-based
strategies, perform on par with each other and outperform the random strategy
by a large margin. Both approaches discover 80% of the top 1% of synergistic
drug pairs in approximately 15 iterations whereas the random strategy
discovers less than 20% in the same number of iterations.

Figure S7B presents the average synergy among queried combinations at each
iteration of the SMO pipeline. While the average synergy using random strategy
is always approximately $10$, model-based strategies query batches for which
the average synergy can be up to 25. After approximately 15 iterations, the
average synergy in the queries starts to decrease, as there are fewer highly
synergistic combinations left to query.

\begin{overpic}[width=195.12767pt]{img/appendix/FigureS7_active_learning_top_zommed_in.png}
\put(2.0,65.0){A} \end{overpic}

\begin{overpic}[width=195.12767pt]{img/appendix/FigureS7_active_learning_average_syn_zoomed_in.png}
\put(45.0,10.0){\includegraphics[clip,width=113.81102pt]{img/appendix/FigureS7_active_learning_average_syn.png}}
\put(2.0,65.0){B} \end{overpic}

Figure S7. Comparison of acquisition functions through in silico SMO
experiments: (A.) Cumulative proportion of the top 1% synergistic combinations
that have been rediscovered by RECOVER, (B.) Average of the true synergies of
the drug pairs that have just been queried at each iteration of the SMO
pipeline, (B, inset) Zoomed out view. Uncertainty estimated via deep
ensembles. Shaded areas correspond to the standard deviation computed over 3
seeds.

###  C.3. Transfer learning: O’Neil 2016 study to NCI-ALMANAC study

An important aspect of the SMO pipeline is its ability to leverage publicly
available data in order to improve performance in a new experimental setting.
To this end, we analyzed the impact of pretraining the model on the O’Neil
database before simulating SMO experiments on a subset of the NCI-ALMANAC
database.

While the out-of-distribution analysis presented in Table 1 showed that
RECOVER does not generalize well to new experimental settings without
adaptation, these experiments demonstrate, quite remarkably, that some latent
knowledge can still be transferred from one experimental setting to another,
resulting in increased performance in the latter setting.

As shown in Figure S8, the model pretrained on O’Neil outperforms the other
model, initialized randomly. The model was first pretrained on the O’Neil
study, thereafter, we simulated the SMO process on the subset of NCI-ALMANAC
consisting of drug pairs for which at least one of the drugs was included in
the O’Neil study. We compared against a model that had not been pretrained.
All models are conditioned on cell line, and two different uncertainty
estimation methods were tested.

In these experiments, we restricted ourselves to cell lines which were covered
in both the O’Neil and NCI-ALMANAC studies. As the MCF7 cell line is not
included in the O’Neil study, we could not restrict ourselves to MCF7 as
usual. Instead, we performed SMO on all overlapping cell lines, resulting in a
space of possible queries which is bigger than in other experiments where only
MCF7 was used. This explains why the rate of discovery is slower in this case.

\begin{overpic}[width=195.12767pt]{img/appendix/FigureS8_effect_of_pretraining_ensemble.png}
\put(2.0,65.0){A} \end{overpic}

\begin{overpic}[width=195.12767pt]{img/appendix/FigureS8_effect_of_pretraining_direct_uncert.png}
\put(2.0,65.0){B} \end{overpic}

Figure S8. Effect of pretraining on the rate of discovery of highly
synergistic combinations via a UCB acquisition function: (A.) Deep Ensemble,
(B.) Direct Uncertainty Estimation. Evaluation performed on the subset of NCI-
ALMANAC consisting of drug pairs for which at least one of the drugs is
included in O’Neil. All overlapping cell lines were included. Shaded areas
correspond to the standard deviation computed over 3 seeds.

###  C.4. Demonstration that an uncertainty driven strategy can outperform a
greedy strategy

When trying to mirror as closely as possible the setting of the in vitro
experiments, we did not notice any significant difference in performance
between the two Model-based acquisition functions Greedy and UCB, as shown in
Figure S7A. However, depending on the model set up, we can demonstrate that
taking uncertainty into account to guide experiments can increase the
performance of the pipeline over a naive Greedy acquisition strategy.

In the following, uncertainty was directly estimated using an DEUP-style
uncertainty predictor, and the task is the prediction of the average Bliss
synergy score (instead of maximum Bliss), which corresponds to the average
over the dose-response matrix of the concentration specific Bliss scores. In
the following, 5 combinations were acquired at a time, instead of 30.

As shown in Figure S9, UCB can outperform Greedy acquisition, demonstrating
the value of taking uncertainty into account in the exploration strategy.

Figure S9. Comparison of acquisition functions through in silico SMO
experiments. Prediction of the average Bliss synergy score. 5 combinations
acquired at a time. Uncertainty estimated using direct uncertainty estimation.
Shaded areas correspond to the standard deviation computed over 3 seeds.