# DEL-Dock: Molecular Docking-Enabled Modeling of DNA-Encoded Libraries

Kirill Shmilovich  [  Benson Chen  [  Theofanis Karaletsos  [
<EMAIL_ADDRESS>Mohammad M. Sultan  [<EMAIL_ADDRESS>
###### Abstract

DNA-Encoded Library (DEL) technology has enabled significant advances in hit
identification by enabling efficient testing of combinatorially-generated
molecular libraries. DEL screens measure protein binding affinity though
sequencing reads of molecules tagged with unique DNA-barcodes that survive a
series of selection experiments. Computational models have been deployed to
learn the latent binding affinities that are correlated to the sequenced count
data; however, this correlation is often obfuscated by various sources of
noise introduced in its complicated data-generation process. In order to
denoise DEL count data and screen for molecules with good binding affinity,
computational models require the correct assumptions in their modeling
structure to capture the correct signals underlying the data. Recent advances
in DEL models have focused on probabilistic formulations of count data, but
existing approaches have thus far been limited to only utilizing 2-D molecule-
level representations. We introduce a new paradigm, DEL-Dock, that combines
ligand-based descriptors with 3-D spatial information from docked protein-
ligand complexes. 3-D spatial information allows our model to learn over the
actual binding modality rather than using only structure-based information of
the ligand. We show that our model is capable of effectively denoising DEL
count data to predict molecule enrichment scores that are better correlated
with experimental binding affinity measurements compared to prior works.
Moreover, by learning over a collection of docked poses we demonstrate that
our model, trained only on DEL data, implicitly learns to perform good docking
pose selection without requiring external supervision from expensive-to-source
protein crystal structures.

University of Chicago] Pritzker School of Molecular Engineering, University of
Chicago, Chicago, Illinois 60637, United States Equal Contribution Insitro]
Insitro, South San Francisco, California 94080, United States Equal
Contribution Insitro] Insitro, South San Francisco, California 94080, United
States Insitro] Insitro, South San Francisco, California 94080, United States

##  1 Introduction

One key component of early drug discovery is hit identification, or finding
molecules with desired activity levels and properties of interest 1. This task
has been commonly approached via high-throughput screening (HTS) techniques,
which involve testing a library of molecules against a biological target of
interest. While traditional screening techniques do not scale well to large
chemical spaces, DNA-encoded library (DEL) technologies provide one avenue
towards the desired scale. For instance, while traditional HTS libraries
typically contain only $\sim$50k-5M compounds, DELs enable screening of
combinatorially large molecular spaces that allow for testing $\sim$1M-5B
compounds in a single tube 2, 3. By capturing a broad chemical diversity
landscape, DEL technology has opened new opportunities in hit identification
4, 5, 6, 7.

DELs are constructed by sequentially assembling molecular building blocks, aka
synthons, into molecules tagged with unique DNA-barcode identifiers (Fig. 1).
Once synthesized, the library is tested for affinity against a protein of
interest through a series of selection experiments. This process, called
panning, typically involves spiking the DEL into a solution of the immobilized
protein, and washing the resulting mixture for multiple rounds. This procedure
leaves members of the library that remain bound to either the target or
matrix, which are subsequently identified using next-generation DNA
sequencing. The resulting data after bioinformatics processing consists of
sparse reads of the DNA and the corresponding molecules; the relative
abundance of the identified molecules is, in theory, a reasonable proxy for
their binding affinities. However, DEL data from panning experiments contain
various sources of noise such as matrix binding, truncation products, unequal
initial loads, replicate and sequencing biases 8, 9, 10, 11. Modeling the data
without considering these significant sources of noise can overfit to spurious
correlations. Consequently, while DELs have proven to be powerful tools in
drug discovery, careful consideration of the appropriate techniques to denoise
DEL data is required to discover reliable signals of the underlying small
molecule binding affinities.

Figure 1: Step-by-step diagram of a DNA-Encoded Library (DEL) panning
experiment.

Several approaches have been developed to tackle DEL modeling from a
computational perspective: these existing methods rely on computing an
enrichment score for each molecule that is representative of its binding
affinity to the protein target. One approach for calculating enrichment scores
involves fitting a Poisson distribution to observed on-target and control
count data for each molecule and then formulating enrichment scores as ratios
of these Poisson parameters 12. A deficiency in this enrichment metric is that
it neglects to include any structural information of the actual DEL molecules.
Since molecular structure represents an important aspect for determining
protein-ligand interactions 13, 14, other methods incorporate different types
of molecular representations into their models. McCloskey et al. 15 use a
Graph Neural Network (GNN) to encode molecules, formulating the learning
problem as multi-class binary classification over the counts. Recent methods
consider hybrid approaches that combine molecular representation learning with
probablistic modeling8, 9, 16. Central to these approaches involves modelling
observed sequencing counts as originating from latent Poisson or Gamma-Poisson
distributions, which are well suited to describe independently distributed
count events. This allows the models to incorporate the inherent noise in the
data-generating process within the model structure. However, these methods
focus on learning only representations of the molecules themselves and omit
the rich representation of the 3-D protein-ligand interactions which are
critical aspects for describing protein binding.

In order to leverage 3-D spatial data, we turn to molecular docking, which has
a long history of being an important tool in structure-based drug discovery
17. A key component of these docking methods is to define a scoring function
that can characterize the binding interactions between a protein and ligand.
To generate candidate poses, these methods use the aforementioned scoring
function to sample ligand conformations with the goal of ultimately inferring
the most likely binding mode. While there have been many advances to docking
models and software, in practice, there are still many pitfalls of these
methods 18. For instance, the difficulty of the docking problem can be
observed through the low empirical correlation between high scoring docked
poses and actual binding affinity 19. Improvements can be realized through
carefully calibrated scoring functions, but tuning these scoring functions is
frequently problem-specific and involves substantial domain knowledge 20, 21,
22. ML-enabled scoring functions trained on databases of crystal structures
have recently demonstrated improved performance in some settings 23, however,
these approaches are ultimately limited by scarce and expensive-to-generate
crystal structures. Even though docking generates noisy ligand conformations,
these docking techniques can provide a useful distribution of likely binding
poses which we will exploit in our new ideas for modeling DEL count data.

Towards the goal of more holistic DEL models, we propose DEL-Dock, a model
that directly learns a joint protein-ligand representation by synthesizing
multi-modal information within a unified probabilistic framework to learn
enrichment scores. Our approach combines molecule-level descriptors with
spatial information from docked protein-ligand complexes to explain sequencing
counts. We delineate contributions from spurious matrix binding and target
protein binding to better separate the signal from noise in the data.
Additionally, our model can learn to rank a collection of poses without
explicit supervision of pose scores by only training on the count data. When
viewed separately, DEL data and docked poses provide noisy signals of binding
affinity, but our DEL-Dock model effectively combines these two modalities to
better extract the learning signals from the data.

In this work, we engage the commonly studied protein human carbonic anhydrase
(CAIX) by training our model on publicly available DEL data generated by Gerry
et al. 12. On a held-out evaluation data set of 3041 molecules with
experimental affinity measurements extracted from the BindingDB 24 web
database, we demonstrate our approach effectively learns enrichment scores
that are well correlated with binding affinity. Moreover, we show that our
model is capable of extracting insights into the determinants of docked
protein-ligand complexes that are most influential for protein binding.

##  2 Methods

Our model combines two different representation modalities, molecule-level
descriptors and docked protein-ligand complexes, to capture the latent aspects
of protein binding from a probabilistic perspective (Fig. 2). The
combinatorial construction of DELs motivates using expressive molecular
representations to capture statistical correlation between the building block
substructures used in DEL synthesis. We use Morgan fingerprint, calculated
with RDKit (version 2020.09.1) 25, as basis for our molecular representations,
which is a standard descriptor for representational problems on small
molecules 26 and also provide the added benefit of their simple construction
and rapid processing. Morgan fingerprints compute a structural bit hash of the
molecule by enumerating $k$-hop substructures about each atom. Since there are
many shared structural features across different molecular compounds, these
fingerprints constitute a simple representation that has demonstrated
remarkable empirical performance throughout cheminformatic domains. We
represent docked protein-ligand poses using a pretrained voxel-based CNN model
from GNINA, which captures spatial relationships by discretizing space into
three-dimentional voxels and leveraging CNNs to learn complex hierarchical
representations 27, 28, 29. The CNN models used in this work was originally
trained on the PDBBind 30 database, capitalizing on this supervised data
source to capture the important features that characterize protein-ligand
interactions.

Figure 2: Schematic illustration of our DEL-Dock neural network architechture
and data flow.

Let $\mathcal{X}$ denote the set of molecules in our data, where each molecule
$x\in\mathcal{X}$ has an associated set of $n$ docked poses
$\\{p_{1},p_{2},...,p_{n}\\}\in\mathcal{P}$ and $c^{\text{matrix}}_{i}\in
C^{\text{matrix}},c^{\text{target}}_{i}\in C^{\text{target}}$ are the $i$th
replicates (repeated experiments) of count data from the beads-only control
and target protein experiments respectively. Additionally, we can define the
following featurization transformations that are used to construct the
molecule and pose embeddings: $\Phi:\mathcal{X}\rightarrow[0,1]^{n_{\phi}}$ is
the function that generates a $n_{\phi}$-bit molecular fingerprint; here, we
use a 2048-bit Morgan fingerprint with radius 3.
$\Psi:\mathcal{X}\times\mathcal{P}\rightarrow\mathbb{R}^{n_{\psi}}$ is the
transformation that outputs an embedding of the molecule and a specific
spatial protein-ligand complex, where we use a pre-trained voxel-based CNN to
perform this transformation.

Let $h^{\text{fps}}=\text{MLP}(\Phi(x))$ be the molecule embedding learned by
our model, which is computed by applying a multilayer perceptron (MLP) to the
fingerprint representation. Individual docked pose embeddings are similarly
computed, with one difference being that we also incorporate the fingerprint
embedding, $h^{\text{pose}}_{i}=\text{MLP}([\Psi(x,p_{i});h^{\text{fps}}])$
into this representation.

To synthesize the set of poses for each molecule we apply a self-attention
layer over the pose embeddings. Following previous work on self-attention and
multiple-instance learning (MIL) 31, we compute attention weights in eq. 1,
where $\big{(}w,W^{U},W^{V}\big{)}$ are learnable weights, $\sigma$ is the
sigmoid activation and $\odot$ is element-wise multiplication. The final
output pose embedding that combines information from the individual input
poses is then computed as a attention-score weighted embedding vector
$h^{\text{pose}}=\frac{1}{n}\sum_{i}a_{i}h^{\text{pose}}_{i}$.

$a_{i}=\frac{\exp\Big{[}w\cdot\big{(}\text{tanh}(W^{U}h^{\text{pose}}_{i})\odot\sigma(W^{V}h^{\text{pose}}_{i})\big{)}\Big{]}}{\sum_{j}\exp\Big{[}w\cdot\big{(}\text{tanh}(W^{U}h^{\text{pose}}_{j})\odot\sigma(W^{V}h^{\text{pose}}_{j})\big{)}\Big{]}}$
(1)

Equipped with these molecule and pose embeddings, our model learns the
contributions of both spurious matrix binding and target protein binding by
predicting latent scores that strive to maximize the likelihood of the
observed data under the model. We make several distinct modeling choices to
mirror our assumptions about the data-generation process that accounts for
various sources of experimental noise.

Matrix Binding is a confounding factor inherent to DEL experiments, since
molecules are prone to binding to the multifarious components comprising the
immobilized matrix in addition to the intended protein target. For each
molecule $x$, we learn a latent matrix binding score
$\lambda^{\text{matrix}}=f(h^{\text{fps}})$. Since matrix binding is not a
function of the protein-ligand pose representation, we enforce that the matrix
binding enrichment remains only a function of the molecule embedding
$h^{\text{fps}}$.

Target Binding is learned through
$\lambda^{\text{target}}=f(h^{\text{pose}})$, jointly utilizing both molecule
and pose representations. This design choice reflects that sequencing counts
from the target protein experiment must be a function of both small molecule
binding to the protein receptor, represented here as featurizations of the
docked protein-ligand complexes, along with promiscuous binding to the
immobilized matrix.

The observed count data for both the control and protein target experiments
can be modeled as originating from underlying Poisson distributions, which
naturally characterize any discrete count data from independently sampled
events. Due to possible sequencing noise we further augment this basic Poisson
model as a zero-inflated probability distribution. This design choice is
motivated by the chance that sparse zero counts in the data could be explained
as an artifact of imperfect sequencing technology, and as a result we directly
incorporate this assumption into the structure of our model. We note that
previous approaches have employed Gamma-Poisson distributions to model DEL
count data in the past 16, but we make the assumption here that replicate data
for the same molecule is sampled from an identical distribution.

$\text{P}(C=c|\lambda,\pi)=\begin{cases}\pi+(1-\pi)e^{-\lambda}&\text{if
}c=0\\\ (1-\pi)\frac{\lambda^{c}e^{-\lambda}}{c!}&\text{if }c>0\\\
\end{cases}$ (2)

$\displaystyle C^{\text{matrix}}$
$\displaystyle\sim\text{ZIP}(\lambda^{\text{matrix}},\pi^{\text{matrix}})$
$\displaystyle\lambda^{\text{matrix}}$
$\displaystyle=\exp(\text{MLP}(h^{\text{fps}}))$ (3) $\displaystyle
C^{\text{target}}$
$\displaystyle\sim\text{ZIP}(\lambda^{\text{matrix}}+\lambda^{\text{target}},\pi^{\text{target}})$
$\displaystyle\lambda^{\text{target}}$
$\displaystyle=\exp(\text{MLP}(h^{\text{pose}}))$ (4)

Let $C$ be distributed as a zero-inflated Poisson (ZIP), with its probability
density function (PDF) defined in eq. 2. Here, $\lambda$ is the rate parameter
of the underlying Poisson distribution, and $\pi$ denotes the occurrence of
choosing the zero distribution, and is taken to be the empirical average (when
$\pi$ = 0, this reduces to the typical Poisson distribution). Empirically, we
estimate these zero-count probabilities from the zero-count frequencies in the
control ($\pi^{\text{matrix}}\approx 0.0075$) and protein target
($\pi^{\text{target}}\approx 0.55$) experiments. Since we model the control
and protein experiments as originating from separate underlying count
distributions, we compute two distinct rate parameters for each ZIP
distribution as shown in eq. 3. The observed target counts is a function of
both matrix binding and binding to the protein target, so the rate parameter
for the target distribution is a function of $\lambda^{\text{matrix}}$ and
$\lambda^{\text{target}}$. We assume an additive functional form for the
latent matrix and target enrichments, which we find works well empirically. We
note however that there could be more motivated methods to parameterize the
interactions of matrix and target binding. The final loss function is then a
typical negative log-likelihood (NLL) loss over the observed counts from both
the control and target experiments eq. 5.

$L=-\sum_{i}\log\big{[}P(c^{\text{matrix}}_{i}|\lambda^{\text{matrix}},\pi^{\text{matrix}})\big{]}-\sum_{j}\log[P(c^{\text{target}}_{j}|\lambda^{\text{target}}+\lambda^{\text{matrix}},\pi^{\text{target}})]$
(5)

###  2.1 DEL Data

To train our model we use publicly available DEL data collected by Gerry et
al. 12. This tri-synthon library consists of $\sim$100k molecules with count
data for panning experiments for the human carbonic anhydrase IX (CAIX)
protein. In addition to on-target counts, the data includes beads-only no-
target controls. Four replicate sets of counts for the protein target
experiments are provided, while two replicates of the control experiments are
provided in this data set. To account for possible noise in different
replicates, we follow previous work and normalize the counts for each target
and control replicate by dividing each count by the sum of counts in that
replicate experiment and then multiplying by 1e6 to re-calibrate the scale of
the counts 16. This data prepossessing provides the interpretation of each
molecule count as a molecular frequency of that molecule within the DEL
library. The processed data set is then used to train our models employing a
80/10/10 - train/validation/test split. Complete training details are further
provided in the Supporting Information.

###  2.2 Evaluation Data

We evaluate the performance of our models on benchmarks using an external set
of affinity measurements of small molecule curated from the BindingDB 24 web
database. We queried binding affinities for the Human Carbonic anhydrase 9
(CAIX) protein target (UniProt: Q16790), and kept only molecules containing
the same atom types as those present in the DEL data set (C, O, N, S, H, I).
This external evaluation data set is composed of 3041 small molecules with
molecular weights ranging from $\sim$25 amu to $\sim$1000 amu and associated
experimental inhibitory constant (Ki) measurements ranging from $\sim$0.15 M
to $\sim$90 pM. We use the median affinity value in the cases where multiple
different affinity measurements were reported for the same molecule. We also
consider a subset of this dataset which consists of the 521 molecules with
molecular weights between 417 amu and 517 amu (Fig. 6 in the Supporting
Information). These molecular weights correspond to the interquartile range
bounding the 10th and 90th percentiles of the molecular weights in training
dataset. This restricted subset presents a more challenging test as
differentiation cannot rely on only extensive properties such as molecular
weight, but must also effectively identify chemical motifs that impact
molecular binding within this tightly bound range of molecular weights. We
notice that simple properties such as molecular weight or benzenesulfonamide
presence, which is known to be an important binding motif for carbonic
anhydrase 12, 32, 33, achieve better baseline performance on the full
evaluation data compared to the restricted subset. These metrics suggest that
this subset is more challenging as predictors must learn beyond these simple
molecular properties to achieve good performance.

###  2.3 Docking

We perform molecular docking to generate a collection of ligand-bound poses to
a target protein of interest for all molecules within in our training and
evaluation data sets. Docking is performed using the GNINA docking software
27, 28, 29 employing the Vina 34 scoring function. All molecules are docked
against CAIX (PDB:5FL4) with the location of the binding pocket determined by
the bound crystal structure ligand (9FK), using the default GNINA settings
defining an $8\times 8\times 8$ Å3 bounding box around this ligand. Initial
three-dimensional conformers for all docked molecules were generated with
RDKit. For each molecule, we obtain 20 docked poses from GNINA using an
exhaustiveness parameter of 50, using the Vina scoring for end-to-end pose
generation. This approach can similarly be performed using AutoDock Vina 34 or
Smina 35 using the Vina scoring function.

##  3 Results

We demonstrate that our model outperforms previous systems on DEL enrichment
prediction by jointly combining topological features from the molecular graph
and the spatial 3-D protein-ligand information, and additionally illustrate
the capability of our model to better rank ligand poses compared to
traditional docking. Our model learns latent binding affinity for each
molecule to both the matrix and the target as the denoised signals compared to
the observed count data. In this interpretation we should expect higher
enrichment scores predicted by our model to be well-correlated with binding
affinity, and therefore provide a useful metric for predicting anticipated
protein binding in virtual screening campaigns.

###  3.1 DEL-Dock outperforms baselines using only docking or molecule-level
descriptors

To evaluate the performance of our model, especially with respect to out-of-
domain protein binding prediction, we first train our model on DEL data
screened against the human carbonic anhydrase IX (CAIX) protein target 12, and
then predict enrichment scores for molecules with externally measured
experimental binding affinities to CAIX. We evaluate performance in this
setting by measuring spearman rank-correlation coefficients between predicted
enrichments and the experimental affinity measurements, which is a metric that
is agnostic to the scale of the values. Our model only restricts the
enrichment scores to be positive quantities, with no specific distributional
constraints, so spearman rank-correlation, which computes a correlation based
only on the ordinal ranking of the predicted enrichments, is well suited for
our test scenario.

Our method, DEL-Dock, which combines information from docked complexes with
molecular descriptors outperforms previous techniques which only utilize one
of these two data modalities (see Table 1). We find that traditional docking
scores alone generated from AutoDock Vina 34 result in the worst overall
correlations, commensurate with previous observations that docking scores
alone are typically not reliable predictors of binding affinity 19.
Performance based on docked poses alone is however greatly improved when re-
scoring the docked poses using pretrained GNINA 36 CNN models. Another set of
baselines we consider are DEL models that rely only on molecular descriptors.
First, we consider a simple model that involves training a random forest (RF)
on the Morgan fingerprints using the enrichment metrics originally formulated
to facilitate analysis of the DEL data set by Gerry et al. 12.

We then similarly predict the enrichment scores of molecules in the held-out
dataset and the correlation with experimental Ki data. This baseline achieves
reasonable performance on both the full and subset evaluation data, especially
given the simplicity of the model. For a more sophisticated baseline, we train
the Graph Neural Network (GNN) model with the DEL-specific loss function
defined by Lim et al. 9 While this approach achieves good performance on the
full evaluation data, correlations on the restricted subset are largely
unchanged for all docking-based and molecular descriptor-based baselines. Our
model which combines docking pose embeddings with molecular fingerprint
representations outperforms all other baselines, with the largest improvements
of $\sim$2$\times$ better spearman correlations than other approaches realized
on the more challenging molecular weight restricted subset. We provide further
ablation studies of our model in the Supporting Information.

Model | Spearman Ki (full) ↓ | Spearman Ki (subset) ↓
---|---|---
Molecular weight | -0.121 | 0.074
Benzenesulfonamide presence | -0.199 | -0.063
Top Vina 34 docking score | -0.068 | 0.119
Top GNINA 36 docking score | -0.279 $\pm$ 0.044 | -0.091 $\pm$ 0.061
| RF trained on enrichment
---
scores from ( Gerry et al. 12)
-0.231 $\pm$ 0.007 | -0.091 $\pm$ 0.012
GNN ( Lim et al. 9) | -0.298 $\pm$ 0.005 | -0.075 $\pm$ 0.011
DEL-Dock (ours) | -0.328 $\pm$ 0.01 | -0.186 $\pm$ 0.01
Table 1: Comparison of spearman rank-correlation coefficients between
predicted affinity scores and experimental inhibition constant ($K_{i}$)
measurements curated from BindingDB 24. Spearman correlations are shown for
the complete 3041-molecule data set (full), and a 521-molecule subset of this
full data set confined to molecular weights between 417-517 amu. This
molecular weight range approximately corresponds to the 10th and 90th
interquartile range of the molecular weights spanned by the DEL data set.
Error bars are reported as standard deviations over five independently
initialized models.

###  3.2 DEL-Dock better ranks molecules with known binding chemical motifs

While our approach displays good prediction accuracy with respect to
experimental binding measurements, we also find our model provides insights
into the structural and chemical factors that influence binding. Compounds
containing benzenesulfonamide have been well established in literature as the
primarily chemical motif that drives small molecule binding to carbonic
anhydrase 12, 32, 33. Though we do not explicitly incorporate this as a
learning signal for our model, we observe that our model is able to learn this
association, visibly predicting sulfonamides within our evaluation data set as
more highly enriched compared molecules which do not contain
benzenesulfonamides (Fig. 3a). Interestingly, we observe a comparatively large
fraction of non-benzenesulfonamides identified as good binders with low
experimental $K_{i}$. The elevated population of highly enriched non-
benzenesulfonamides in this data set could be an artifact of bias in
scientific literature. Our model is ultimately trained on DEL data and
therefore is expected to reflect underlying biases and idiosyncrasies of the
data generation process. The most notable difference lies in that DEL
experiments are only capable of measuring on-DNA binding, while the evaluation
data are measurements of off-DNA binding. Nevertheless, the clear delineation
of benzenesulfonamides in our predicted enrichments provides good post-hoc
evidence that our model correctly identifies this important binding motif for
this protein target.

Figure 3: Analysis of DEL-Dock model predictions on our evaluation data set
composed of experimental affinity (Ki) measurements. (a) Parity plot of our
model predicted enrichments and ground truth Ki measurements delineated by
benzenesulfonamide presence, with 1581 benzenesulfonamide-containing molecules
and 1460 non-benzenesulfonamide. (b) Distribution of zinc-sulfonamide
distances for the top-ranked pose of each benzenesulfonamide in our evaluation
data set identified by the AutoDock Vina scoring function 34, GNINA pose
selection 36, and our model predicted attention scores. (c) Comparison in the
fraction of top ranked poses identified with with zinc-sulfonamide distance
below a threshold distance between $\sim$2-12Å. This can be interpreted as the
CDF of the distributions in (b) with the appropriate normalization.

An important structural component of benzenesulfonamides binding to carbonic
anhydrase is coordination of the sulfonamide group with the zinc ion buried
within the active site 12, 32, 33. In the vast majority of cases, one would
then expect docking scoring functions to highly score poses that reflect this
anticipated binding mode. As our model performs self-attention over pose
embeddings, which are used to learn molecules’ enrichment scores, we can
interpret the magnitude of the attention probabilities as the importance
weight of that particular pose. Shown in Fig. 3b is the distribution of zinc-
sulfonamide distances for the top-selected docked pose comparing AutoDock
Vina, GNINA, and our method for all 1581 benzenesulfonamides-containing
molecules in our evaluation data set. An alternate view of this data is
presented as the fraction of top-selected poses with zinc-sulfonamide
distances below a distance threshold (Fig. 3c), which can effectively be
interpreted as the cumulative distribution function (CDF) of the appropriately
normalized associated probability distribution function (PDF) in Fig. 3b.

The AutoDock Vina scoring function exhibits largest spread of zinc-sulfonamide
distances, and as a result identifying a comparatively large fraction of poses
as incorrectly coordinated. GNINA pose selection performs significantly better
in this setting, identifying a larger fraction of well-coordinated poses with
low zinc-sulfonamide distance. We find our method ultimately correctly
coordinates the largest proportion of poses when compared to AutoDock Vina or
GNINA. We note that our approach for binding pose selection is markedly
different than the approach taken by GNINA, which involves a separate pose
scoring head trained to identify poses with low-RMSD to ground truth crystal
structures. Our attention scores on the other hand are effectively latent
variables trained only via the auxiliary task of modeling DEL data. The
benefit of our approach is that we can learn to identify good poses in an
unsupervised manner, without requiring scarce and expensive crystal structures
to serve as the source of supervision for pose selection.

###  3.3 DEL-Dock offers interpretability through its attention mechanisms

Lastly, we further demonstrate the interpretability of or model by examining
the distribution of attention scores learned by our model for a specific
molecule (Fig. 4). For this molecule, only 7 out of 20 docked poses correctly
coordinate the sulfonamide group with the zinc ion buried in the protein
active site. Our model appropriately identifies this binding mode and learns
attention scores that more favorably rank these 7 correctly coordinated poses
(Fig. 4). The top-three ranked poses (Fig. 4a) by our model have very similar
conformations, each exhibiting zinc-sulfonamide coordination, and differing
only in the orientation of the terminal benzene ring that is distant from the
active site. The other poses that show zinc-sulfonamide coordination (Fig.
4b-d) are also ranked highly by our model, however, these poses exhibit less
favorable conformations in several ways. For instance, the conformation in
Fig. 4b is more exposed, and less protected by the protein. Finally, our model
in general more poorly ranks poses that display incorrect zinc-sulfonamide
coordination (Fig. 4e). These conformations typically have the terminal
benzene ring inserted into the active site. Also, these poses reveal why zinc-
sulfonamide distances alone can be a deceiving metric as some poses are
capable of achieving low zinc-sulfonamide distances ($\sim$3 Å) due to the
molecule “curling in” on itself within the active site. Nevertheless, our
model recognizes this spurious binding mode and poorly ranks these poses with
comparatively low attention scores, even though AutoDock Vina highly ranks
many of these bad poses. Overall, we find the hierarchy of pose rankings by
our model to be commensurate with anticipated binding behavior for this
protein target.

Figure 4: Analysis of pose attention scores for a representative molecule in
our evaluation data set. (left) Our model predicted pose attention scores
plotted against the zinc-sulfonamide distance of the docked pose and colored
according to the ranking determined by the AutoDock Vina scoring function.
(right) Different protein-ligand complexes are visualized to show that our
model highly ranks the conformers with zinc-sulfonamide coordination (a-d),
while the conformers without the correct coordination are ranked lower (e).

##  4 Conclusions

In this work we present an approach for modeling DEL data that combines
docking-based and molecular descriptor-based data modalities. Our approach
involves predicting two interleaved quantities, enrichment scores, that
explain the sequencing counts of the panning experiment measurements for both
the on-target protein and the off-target control beads. We evaluate our method
by first training on DEL data screened against the human carbonic anhydrase
(CAIX) 12 protein target, and then predicting binding for unseen molecules
with external experimental constant of inhibition ($K_{i}$) affinity
measurements curated from the BindingDB 24 web database. For this prediction
task we find our approach outperforms previous docking and DEL modeling
techniques that only use either docked poses or molecular descriptor
information alone. Furthermore, a critical component of our model involves
performing self-attention over pose embedding, in order to learn over the set
of possible poses. Analyzing these latent attention scores, we find our model
effectively identifies good docked poses. Compared to docking pose selection
using either AutoDock Vina 34 or GNINA 27, 28, 29, 36, our model more reliably
selects poses displaying the appropriate zinc-sulfonamide coordination–which
is known to be the predominant binding mode for carbonic anhydrase 12, 32, 33.
Our model is interestingly capable of learning good pose selection in an
unsupervised manner, training only on the voluminous DEL data rather than
requiring crystal structures to serve as the source of supervision.

There are some limitations to our approach, however. Our approach assumes that
we can obtain reasonably accurate protein-ligand docking poses, which is not
always true. We can only utilize this model for proteins which have 3-D
crystal structures; and even if protein crystal structures are available, they
may not always be accurate. Additionally, we are limited by the quality of the
docking software, which may fail to capture the actual correct binding modes,
or ranking the good binding modes highly.

Our work focuses on introducing the concept of utilizing docked poses to
improve DEL models, but there are several avenues for future work. While we
use Morgan fingerprints as our molecule featurizer, deep learning approaches
have the potential to generate more expressive representations, such as Graph
Neural Networks (GNNs) that have demonstrated excellent predictive performance
on many molecular property prediction tasks 37. In our application,
equivariant GNNs could also be used as an alternative featurizer to CNNs for
embedding docked poses, which would provide the added benefit of producing
explicitly roto-translationally symmetric representations 38. We only select
the top ranking poses from docking, but providing a more diverse set of poses
could be more useful instead. Future research directions could also leverage
our approach to unsupervised pose selection for downstream free energy
calculations in larger-scale virtual screening campaigns. Our approach paves
the way for multi-modality modeling of DEL data and unsupervised learning of
bespoke, DEL-conditioned, scoring functions.

##  5 Acknowledgements

We would like to thank Nathaniel Stanley and Sam Mun at Insitro for helpful
discussions and advice related to Docking software and usage. We would also
like to thank Daphne Koller, Robert Hilgraf, Nathaniel Stanley, Fiorella
Ruggiu and Yujia Bao at Insitro for providing general feedback and review of
our work. Lastly, we would like to thank Patrick Conrad at Insitro for helping
us with the public code release.

##  6 Data Availability

We use publicly available data from Gerry et al. 12, and our code is publicly
available at: https://github.com/insitro/insitro-research.

##  7 Supporting Information

###  7.1 Training settings

Featurizations for the docked poses are generated using pre-trained GNINA
models provided in gnina-torch 28, 29. The dense variant of the GNINA models
composed of densely connected 3D residual CNN blocks introduced in Francoeur
et al. 36 are used to generate 224-dimentional embeddings of each docked pose.
Morgan fingerprints 26 for each molecule are calculated using RDKit with a
radius of 3 embedded into a 2048 dimensional bit-vector.

All models are trained end-to-end using mini-batch gradient decent with the
Adam optimizer 39 and coefficients for the running averages of
$\beta_{1}=0.95$ and $\beta_{2}=0.999$. A batch size of 64 is used with an
initial learning rate of $1\times 10^{-4}$ and a linearly decaying learning
rate scheduler where the learning rate is decayed by a factor of
$\gamma^{\frac{1}{n_{steps}}}$ every batch. For our learning rate scheduler we
use $\gamma=0.1$ and $n_{steps}=1250$, which corresponds to a $10\times$
reduction in the learning rate after $1250$ batches. We also apply gradient
clipping, where gradient norms are clipped to a maximum value of 0.1. During
training we maintain an exponential moving average over our model parameters
which are updated each step with a decay rate of 0.999. This exponential
moving average version of the model parameters is then used for evaluation and
throughout all inference tasks. Throughout our model we use LeakyReLU
activation functions with a negative slope constant of $1\times 10^{-2}$,
except for the final activation function applied to the output logits
corresponding to the matrix and target enrichment scores where we apply an
exponential function as our terminal activation. A hidden dimensionality of
256 is used within MLP layers in our network. The residual MLP layers, which
are responsible for processing the Morgan fingerprints along with the CNN
features and embeddings (Fig. 2 in the main text), are composed of 2
residually connected MLP layers using dropout with a probability of 0.5. Our
model in sum is composed of $\sim$1M parameters and is trained for 8 epochs on
a single NVIDIA T4 GPU. A PyTorch implementation of our model that makes use
of PyTorch Lightning 40 and Pyro 41 is publicly available at:
https://github.com/insitro/insitro-research.

###  7.2 Ablations

We present here a number of ablations on our model, exploring some different
architectural components and design choices (Table 2). First, we see that
training models with only fingerprint representations, without incorporating
any information from the docked poses, results in a marked decrease in
performance. On the other hand models, trained using only CNN representations
perform much better and display comparable performance to only GNINA pre-
trained models (Table 1). This represents an intuitive result, as this
training setting is effectively equivalent to fine-tuning GNINA using a multi-
instance learning over the pose representations. Interestingly, we find that
training on the CNN features alone already achieves good binding pose
selection based on the latent attention scores, a feature our model is
evidently capable of learning in isolation of the fingerprint representations.
Also shown as a baseline is the MLP network trained using the bespoke loss
function for modelling DEL data presented by Lim et al. 9. This approach
represents the lower performing of the two architectures explored by Lim et
al. 9, the other being the GNN architecture presented as a baseline in the
main text (Table 1).

Model | Spearman Ki (full) ↓ | Spearman Ki (subset) ↓
---|---|---
Only fingerprints | -0.191 $\pm$ 0.005 | -0.083 $\pm$ 0.019
Only CNN | -0.287 $\pm$ 0.005 | -0.124 $\pm$ 0.006
MLP from Lim et al. 9 | -0.244 $\pm$ 0.004 | -0.076 $\pm$ 0.017
| Without zero-inflated
---
distribution
-0.26 $\pm$ 0.02 | -0.08 $\pm$ 0.03
| End-to-end voxels
---
with frozen CNN
-0.278 $\pm$ 0.022 | -0.16 $\pm$ 0.03
Table 2: Model ablations and other baselines. Error bars are calculated as
standard deviations over five independently initialized models.

Interestingly using a zero-inflated loss appears to be critical to our
performance, resulting in $\sim$25% increase in spearman rank-correlation on
the full evaluation set and greater than a 2$\times$ increase on the subset.
We suspect this performance jump could be related to the disparity in zero-
counts between the control and on-target experiment: the control experiments
have a zero-count frequency of $\sim$0.75% while the protein target
experiments have a zero-count frequency of $\sim$55%. Using a zero-inflated
distribution could provide our model more flexibility to explain zero-counts
as an artifact of the data generation process, rather than an outcome of poor
protein binding.

Instead of using pre-computed CNN features from GNINA we also explored
training our model directly from the voxel representations using frozen CNN
featurizers. The benefit of this approach is the ability to use data
augmentation via radom rotations and translations to implicitly enforce that
the learned CNN embeddings remain roto-translationally equivariant. While we
notice performance on the evaluation subset is comparable with our trained on
pre-computed CNN features (Table 1), the performance on the full data set is
slightly reduced. We suspect this result could be due to the computational
challenges of using voxelized representations. In particular, when training
over many docked poses (in our case 20 poses per molecule) our batch size is
effectively 20$\times$ larger – which presents a significant memory bottleneck
as the voxel representation requires storing a
48$\times$48$\times$48$\times$28 molecular grids (three dimensions
discretizing space, and one for different atom types). Furthermore, our pre-
computed features are already being generated with pre-trained CNN featurizers
that have been trained using data augmentation, albeit on PDBBind for the
separate task of affinity and pose prediction. Nevertheless, we certainly
expect improved performances could still be achieved for these full
differentiable training approaches given the appropriate compute resources and
further tuning.

Lastly, presented in Table 3 is a comparison of spearman rank-correlation
performances training on variable numbers of poses. For each model the top-$k$
poses generated via docking are used for training. Performance tends to
generally improve with increasing number of poses used for training, with the
largest difference in improvements realized on the molecular weight restricted
subset. Beyond $\sim$10 poses appears to result in diminishing returns, in
comparison to the jumpy in improvements seen from 2 $\rightarrow$ 10 poses.

Number of training poses | Spearman Ki (full) ↓ | Spearman Ki (subset) ↓
---|---|---
2 poses | -0.278 $\pm$ 0.011 | -0.112 $\pm$ 0.023
5 poses | -0.304 $\pm$ 0.01 | -0.15 $\pm$ 0.02
10 poses | -0.318 $\pm$ 0.007 | -0.175 $\pm$ 0.02
15 poses | -0.324 $\pm$ 0.008 | -0.182 $\pm$ 0.014
20 poses | -0.328 $\pm$ 0.009 | -0.186 $\pm$ 0.013
Table 3: Model ablations training on different numbers of docked poses. Error
bars are calculated as standard deviations over five independently initialized
models.

###  7.3 Supplementary Figures

Shown in Fig. 5 is a TSNE embedding of the DEL data set alongside our
evaluation data. This TSNE embedding is generated by representing each
molecule with a concatenation of three fingerprint representations: a
2048-dimensional Morgan fingerprint with a radius of 3, 167-dimensional MACCS
(Molecular ACCess System) fingerprint, and finally a 2048-dimensional atom
pair fingerprint. All fingerprints are calculated using RDKit. Sci-kit learn
is then used to generate the TSNE embedding using a tanimoto similarity metric
with a perplexity of 30 trained on the combined DEL and evaluation data. We
notice the evaluation data is largely isolated from the DEL data in this TSNE
embedding, serving as an indication that our evaluation data is markedly
different, or out of domain, than the DEL data used in training our models.

Figure 5: TSNE embedding our of DEL and evaluation data. Molecular
representations for this TSNE embedding are generated as a concatenation of
three fingerprint representations: Morgan fingerprints, MACCS fingerprints,
and atom pair fingerprints. Figure 6: Comparison of the distribution of
molecular weights between the DEL data set and (a) the full evaluation data
set and (b) the 417-517 amu subset of the evaluation data set. Distributions
are generated as a Kernel Density Estimate (KDE) plot as implemented in
seaborn 42.

Shown in Fig. 7 are the distributions of zinc-sulfonamide distances throughout
the top-five ranking poses as identified by our DEL-DOCK model attention
scores, GNINA pose selection, and the AutoDock Vina scoring function. We
notice that the highly ranked poses by our model attribute more density in the
closely separated regime under $\sim$4Å than GNINA or Vina, and as a direct
consequence of this we see fewer poses selected by our model showing large
separations between $\sim$4Å - $\sim$13Å.

Figure 7: Distribution of Zinc-Sulfonamide distances throughout the top-five
ranked poses by our DEL-DOCK model attention scores, GNINA pose selection
score, and AutoDock Vina scoring function.

## References

  * Hughes et al. 2011 Hughes, J. P.; Rees, S.; Kalindjian, S. B.; Philpott, K. L. Principles of early drug discovery. _British journal of pharmacology_ 2011, _162_ , 1239–1249
  * Satz et al. 2022 Satz, A. L.; Brunschweiger, A.; Flanagan, M. E.; Gloger, A.; Hansen, N. J.; Kuai, L.; Kunig, V. B.; Lu, X.; Madsen, D.; Marcaurelle, L. A., et al. DNA-encoded chemical libraries. _Nature Reviews Methods Primers_ 2022, _2_ , 1–17
  * Sunkari et al. 2021 Sunkari, Y. K.; Siripuram, V. K.; Nguyen, T.-L.; Flajolet, M. High-power screening (HPS) empowered by DNA-encoded libraries. _Trends in Pharmacological Sciences_ 2021,
  * Clark et al. 2009 Clark, M. A.; Acharya, R. A.; Arico-Muendel, C. C.; Belyanskaya, S. L.; Benjamin, D. R.; Carlson, N. R.; Centrella, P. A.; Chiu, C. H.; Creaser, S. P.; Cuozzo, J. W., et al. Design, synthesis and selection of DNA-encoded small-molecule libraries. _Nature chemical biology_ 2009, _5_ , 647–654
  * Kleiner et al. 2011 Kleiner, R. E.; Dumelin, C. E.; Liu, D. R. Small-molecule discovery from DNA-encoded chemical libraries. _Chemical Society Reviews_ 2011, _40_ , 5707–5717
  * Goodnow et al. 2017 Goodnow, R. A.; Dumelin, C. E.; Keefe, A. D. DNA-encoded chemistry: enabling the deeper sampling of chemical space. _Nature Reviews Drug Discovery_ 2017, _16_ , 131–147
  * Flood et al. 2020 Flood, D. T.; Kingston, C.; Vantourout, J. C.; Dawson, P. E.; Baran, P. S. DNA encoded libraries: a visitor’s guide. _Israel Journal of Chemistry_ 2020, _60_ , 268–280
  * Binder et al. 2022 Binder, P.; Lawler, M.; Grady, L.; Carlson, N.; Leelananda, S.; Belyanskaya, S.; Franklin, J.; Tilmans, N.; Palacci, H. Partial Product Aware Machine Learning on DNA-Encoded Libraries. _arXiv preprint arXiv:2205.08020_ 2022,
  * Lim et al. 2022 Lim, K. S.; Reidenbach, A. G.; Hua, B. K.; Mason, J. W.; Gerry, C. J.; Clemons, P. A.; Coley, C. W. Machine learning on DNA-encoded library count data using an uncertainty-aware probabilistic loss function. _Journal of Chemical Information and Modeling_ 2022,
  * Zhu et al. 2021 Zhu, H.; Foley, T. L.; Montgomery, J. I.; Stanton, R. V. Understanding Data Noise and Uncertainty through Analysis of Replicate Samples in DNA-Encoded Library Selection. _Journal of Chemical Information and Modeling_ 2021, _62_ , 2239–2247
  * Kómár and Kalinic 2020 Kómár, P.; Kalinic, M. Denoising DNA encoded library screens with sparse learning. _ACS Combinatorial Science_ 2020, _22_ , 410–421
  * Gerry et al. 2019 Gerry, C. J.; Wawer, M. J.; Clemons, P. A.; Schreiber, S. L. DNA barcoding a complete matrix of stereoisomeric small molecules. _Journal of the American Chemical Society_ 2019, _141_ , 10225–10235
  * Jones et al. 2021 Jones, D.; Kim, H.; Zhang, X.; Zemla, A.; Stevenson, G.; Bennett, W. D.; Kirshner, D.; Wong, S. E.; Lightstone, F. C.; Allen, J. E. Improved protein–ligand binding affinity prediction with structure-based deep fusion inference. _Journal of chemical information and modeling_ 2021, _61_ , 1583–1592
  * Hwang et al. 2017 Hwang, H.; Dey, F.; Petrey, D.; Honig, B. Structure-based prediction of ligand–protein interactions on a genome-wide scale. _Proceedings of the National Academy of Sciences_ 2017, _114_ , 13685–13690
  * McCloskey et al. 2020 McCloskey, K.; Sigel, E. A.; Kearnes, S.; Xue, L.; Tian, X.; Moccia, D.; Gikunju, D.; Bazzaz, S.; Chan, B.; Clark, M. A., et al. Machine learning on DNA-encoded libraries: a new paradigm for hit finding. _Journal of Medicinal Chemistry_ 2020, _63_ , 8857–8866
  * Ma et al. 2021 Ma, R.; Dreiman, G. H.; Ruggiu, F.; Riesselman, A. J.; Liu, B.; James, K.; Sultan, M.; Koller, D. Regression modeling on DNA encoded libraries. NeurIPS 2021 AI for Science Workshop. 2021
  * Meng et al. 2011 Meng, X.-Y.; Zhang, H.-X.; Mezei, M.; Cui, M. Molecular docking: a powerful approach for structure-based drug discovery. _Current computer-aided drug design_ 2011, _7_ , 146–157
  * Chen 2015 Chen, Y.-C. Beware of docking! _Trends in pharmacological sciences_ 2015, _36_ , 78–95
  * Gupta et al. 2018 Gupta, M.; Sharma, R.; Kumar, A. Docking techniques in pharmacology: How much promising? _Computational biology and chemistry_ 2018, _76_ , 210–217
  * Li et al. 2019 Li, J.; Fu, A.; Zhang, L. An overview of scoring functions used for protein–ligand interactions in molecular docking. _Interdisciplinary Sciences: Computational Life Sciences_ 2019, _11_ , 320–328
  * Pham and Jain 2008 Pham, T. A.; Jain, A. N. Customizing scoring functions for docking. _Journal of Computer-Aided Molecular Design_ 2008, _22_ , 269–286
  * Jain 1996 Jain, A. N. Scoring noncovalent protein-ligand interactions: a continuous differentiable function tuned to compute binding affinities. _Journal of computer-aided molecular design_ 1996, _10_ , 427–440
  * Wallach et al. 2015 Wallach, I.; Dzamba, M.; Heifets, A. AtomNet: a deep convolutional neural network for bioactivity prediction in structure-based drug discovery. _arXiv preprint arXiv:1510.02855_ 2015,
  * Liu et al. 2007 Liu, T.; Lin, Y.; Wen, X.; Jorissen, R. N.; Gilson, M. K. BindingDB: a web-accessible database of experimentally determined protein–ligand binding affinities. _Nucleic acids research_ 2007, _35_ , D198–D201
  * Landrum et al. 2020 Landrum, G. et al. rdkit/rdkit: 2020_09_1 (Q3 2020) Release. 2020; https://doi.org/10.5281/zenodo.4107869
  * Rogers and Hahn 2010 Rogers, D.; Hahn, M. Extended-connectivity fingerprints. _Journal of chemical information and modeling_ 2010, _50_ , 742–754
  * McNutt et al. 2021 McNutt, A. T.; Francoeur, P.; Aggarwal, R.; Masuda, T.; Meli, R.; Ragoza, M.; Sunseri, J.; Koes, D. R. GNINA 1.0: molecular docking with deep learning. _Journal of cheminformatics_ 2021, _13_ , 1–20
  * Ragoza et al. 2017 Ragoza, M.; Hochuli, J.; Idrobo, E.; Sunseri, J.; Koes, D. R. Protein–ligand scoring with convolutional neural networks. _Journal of chemical information and modeling_ 2017, _57_ , 942–957
  * Sunseri and Koes 2020 Sunseri, J.; Koes, D. R. Libmolgrid: graphics processing unit accelerated molecular gridding for deep learning applications. _Journal of chemical information and modeling_ 2020, _60_ , 1079–1084
  * Liu et al. 2017 Liu, Z.; Su, M.; Han, L.; Liu, J.; Yang, Q.; Li, Y.; Wang, R. Forging the basis for developing protein–ligand interaction scoring functions. _Accounts of chemical research_ 2017, _50_ , 302–309
  * Ilse et al. 2018 Ilse, M.; Tomczak, J.; Welling, M. Attention-based deep multiple instance learning. International conference on machine learning. 2018; pp 2127–2136
  * Li et al. 2018 Li, Y.; De Luca, R.; Cazzamalli, S.; Pretto, F.; Bajic, D.; Scheuermann, J.; Neri, D. Versatile protein recognition by the encoded display of multiple chemical elements on a constant macrocyclic scaffold. _Nature chemistry_ 2018, _10_ , 441–448
  * Buller et al. 2011 Buller, F.; Steiner, M.; Frey, K.; Mircsof, D.; Scheuermann, J.; Kalisch, M.; Bühlmann, P.; Supuran, C. T.; Neri, D. Selection of carbonic anhydrase IX inhibitors from one million DNA-encoded compounds. _ACS chemical biology_ 2011, _6_ , 336–344
  * Trott and Olson 2010 Trott, O.; Olson, A. J. AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. _Journal of computational chemistry_ 2010, _31_ , 455–461
  * Koes et al. 2013 Koes, D. R.; Baumgartner, M. P.; Camacho, C. J. Lessons learned in empirical scoring with smina from the CSAR 2011 benchmarking exercise. _Journal of chemical information and modeling_ 2013, _53_ , 1893–1904
  * Francoeur et al. 2020 Francoeur, P. G.; Masuda, T.; Sunseri, J.; Jia, A.; Iovanisci, R. B.; Snyder, I.; Koes, D. R. Three-dimensional convolutional neural networks and a cross-docked data set for structure-based drug design. _Journal of chemical information and modeling_ 2020, _60_ , 4200–4215
  * Stokes et al. 2020 Stokes, J. M.; Yang, K.; Swanson, K.; Jin, W.; Cubillos-Ruiz, A.; Donghia, N. M.; MacNair, C. R.; French, S.; Carfrae, L. A.; Bloom-Ackermann, Z., et al. A deep learning approach to antibiotic discovery. _Cell_ 2020, _180_ , 688–702
  * Stärk et al. 2022 Stärk, H.; Ganea, O.; Pattanaik, L.; Barzilay, R.; Jaakkola, T. Equibind: Geometric deep learning for drug binding structure prediction. International Conference on Machine Learning. 2022; pp 20503–20521
  * Kingma and Ba 2014 Kingma, D. P.; Ba, J. Adam: A method for stochastic optimization. _arXiv preprint arXiv:1412.6980_ 2014,
  * Falcon et al. 2020 Falcon, W. et al. PyTorchLightning/pytorch-lightning: 0.7.6 release. 2020; https://doi.org/10.5281/zenodo.3828935
  * Bingham et al. 2019 Bingham, E.; Chen, J. P.; Jankowiak, M.; Obermeyer, F.; Pradhan, N.; Karaletsos, T.; Singh, R.; Szerlip, P.; Horsfall, P.; Goodman, N. D. Pyro: Deep universal probabilistic programming. _The Journal of Machine Learning Research_ 2019, _20_ , 973–978
  * Waskom 2021 Waskom, M. L. seaborn: statistical data visualization. _Journal of Open Source Software_ 2021, _6_ , 3021