DiffDock: Diffusion Steps, Twists, and Turns for Molecular Docking
Authors: Gabriele Corso
Hannes Stärk
Bowen Jing
Regina Barzilay
Tommi S. Jaakkola
Venue: ICLR 2023 poster (2023)
Keywords: molecular docking
protein-ligand binding
diffusion models
score-based models
molecular structure
equivariance
geometric deep learning

Abstract: Predicting the binding structure of a small molecule ligand to a protein---a task known as molecular docking---is critical to drug design. Recent deep learning methods that treat docking as a regression problem have decreased runtime compared to traditional search-based methods but have yet to offer substantial improvements in accuracy. We instead frame molecular docking as a generative modeling problem and develop DiffDock, a diffusion generative model over the non-Euclidean manifold of ligand poses. To do so, we map this manifold to the product space of the degrees of freedom (translational, rotational, and torsional) involved in docking and develop an efficient diffusion process on this space. Empirically, DiffDock obtains a 38% top-1 success rate (RMSD<2A) on PDBBind, significantly outperforming the previous state-of-the-art of traditional docking (23%) and deep learning (20%) methods. Moreover, while previous methods are not able to dock on computationally folded structures (maximum accuracy 10.4%), DiffDock maintains significantly higher precision (21.7%). Finally, DiffDock has fast inference times and provides confidence estimates with high selective accuracy. 

Published as a conference paper at ICLR 2023
DIFFDOCK: DIFFUSION STEPS, TWISTS, AND TURNS FOR MOLECULAR DOCKING
Gabriele Corso", Hannes Stärk*, Bowen Jing', Regina Barzilay & Tommi Jaakkola CSAIL,, Massachusetts Institute of Technology
ABSTRACT
Predicting the binding structure of a small molecule ligand to a protein-a task known as molecular docking- critical to drug design. Recent deep learning methods that treat docking regression problem have decreased runtime com- pared to traditional search-based methods but have yet to offer substantial im- provements in accuracy. We instead frame molecular docking modeling problem and develop DIFFDOCK, a diffusion generative model over the non-Euclidean manifold of ligand poses. To do so, we map this manifold the product space of the degrees of freedom (translational, rotational, and torsional) involved in docking and develop an efficient diffusion process on this space. Em- pirically, DIFFDOCK obtains a 38% top-1 success rate (RMSD<2Á) on PDB- Bind, significantly outperforming the previous atate-of-the-art of traditional dock- ing (23%) and deep learning (20%) methods. Moreover, while previous methods are not able to dock on computationally folded structures (maximum accuracy 10.4%), DurFDOCK maintains significantly higher precision (21.7%). Finally, DIFFDOCK has fast inference times and provides confidence estimates with high selective accuracy.
generative
INTRODUCTION
The biological functions of proteins can be modulated by small molecule ligands (such as drugs) binding to them. Thus, a crucial task in computational drug design molecular Brckkin-preedictng the position, orientation, and conformation of ligand when bound to a target protein-from which the effect of the ligand (if any) might be inferred. Traditional approaches for docking [Trott & Olson, 2010; Halgren et al., 2004] rely on scoring-functions functions that estimate the correctness a proposed structure or pose, and an optimization algorithm that searches for the global maximum of the scoring function. However, since the search space is vast and the landscape of the scoring functions rugged, these methods tend to be too slow and inaccurate, especially for tigh-throughhut workflows.
Recent works [Stärk et al., 2022; Lu et al., 2022] have developed deep learning models to predict the binding pose in one shot, treating docking as a regression problem. While these methods are much faster than traditional search- based methods, they have yet to demonstrate significant improvements in accuracy. We argue that this may because the regression based paradigm corresponds imper- fectly with the objectives of molecular docking, which is reflected in the fact that standard accuracy metries resemble the likelihood of the data under the predictive model rather than a regression loss. We thus frame molecular docking as generative modeling problem- given a ligand and target protein structure, we Iearn a distribution over ligand poses.
To this end, develop DIFFDOCK, diffusion generative model (DGM) over the space of ligand poses for molecular docking. We define a diffusion process over the degrees freedom involved in docking: the position of the ligand relative to the protein (locating the binding pocket), its orientation in the pocket, and the torsion angles describing its conformation. DIFFDOCK samples poses by run- ning the learned (reverse) diffusion process, which iteratively transforms uninformed, noisy prior distribution over ligand poses into the learned model distribution (Figure 1). Intuitively, this pro- cess can be viewed the progressive refinement of random poses via updates of their translations, rotations, and torsion angles.
"Equal contribution. Correspondence {gcorso, hatark, bjing]ami t edu.

---

Published as a conference paper at ICLR 2023
ligand protein
ranked poses confidence score
DIFFDOCK
reverse diffusion over translations, rotations and torsions
t=T
t=0
2
Figure 1: Overview of DIFFDOCK. Left: The model takes as input the separate ligand and protein structures. Center: Randomly sampled initial poses are denoised via reverse diffusion over trans- lational, rotational, and torsional degrees freedom. Right: The sampled poses are ranked by the confidence model to produce final prediction and confidence score.
While DGMs have been applied to other problems in molecular machine learning [Xu et al., 2021; Jing et al., 2022; Hoogeboom et al., 2022], existing approaches are ill-suited for molecular docking, where the space of ligand poses is (m -)dimennion submanifold M R3",, where and m are, respectively, the number of atoms and torsion angles. To develop DIFFDOCK, we recognize that the docking degrees of freedom define M as the space of poses accessible via set of allowed ligand pose transformations. We use this idea to map elements M to the product space of the groups corresponding to those transformations, where DGM can be developed and trained efficiently.
As applications of docking models often require only fixed number of predictions and a confidence score over these, train confidence model provide confidence estimates for the poses sampled from the DGM and to pick out the most likely sample. This two-step process can be viewed as an intermediate approach between brute-force search and one-shot prediction: we retain the ability to consider and compare multiple poses without incurring the difficulties of high-dimensional search.
Empirically, on the standard blind docking benchmark PDBBind, DIFFDOCK achieves 38% of top- predictions with ligand root mean square distance (RMSD) below 2', nearly doubling the perfor- mance the previous state-of-the-art deep learning model (20%). DIFFDOCK significantly out- performs even state-of-the-arr search-based methods (23%), while still being 3 to 12 times faster on GPU, Moreover, it provides an accurate confidence score of its predictions, obtaining 83% RMSD<2A its most confident third of the previously unseen complexes.
We further evaluate the methods on structures generated by ESMFold [Lin et al., 2022]. Our results confirm previous analyses [Wong et al., 2022] that showed that existing methods are not capable of docking against these approximate apo-structures (RMSD<2Á equal or below 10%). Instead, without further training, DIFFDOCK places 22% of its top-1 predictions within 2 opening the way for the revolution brought by accurate protein folding methods in the modeling of protein-ligand interactions.
To summarize, the main contributions of this work are:
1. We frame the molecular docking task as a generative problem and highlight the issues with previous deep learning approaches.
2. We formulate a novel diffusion process over ligand poses corresponding to the degrees of freedom involved in molecular docking.
3. We achieve a new state-of-the-art 38% top-1 prediction with RMSD<2A on PDBBind blind docking benchmark, considerably surpassing the previous best search-based based (23%) and deep learning methods (20%).
4. Using ESMFold to generate approximate protein apo-structures, we show that our method places top-1 prediction with RMSD<2Á on 28% of the complexes, nearly tripling the accuracy of the most accurate baseline.
2

---

Published as a conference paper at ICLR 2023
2 BACKGROUND AND RELATED WORK
Molecular docking. The molecular docking task is usually divided between known-pocket and blind docking. Known-pocket docking algorithms receive as input the position on the protein where the molecule will bind (the binding pocket) and only have to find the correct orientation and confor- mation. Blind docking instead does not assume any prior knowledge about the binding pocket: in this work, we will focus this general setting. Docking methods typically assume the knowledge of the protein holo-structure (bound), this assumption is in many real-world applications unrealis- tic, therefore, we evaluate methods both with holo and computationally generated apo-structures (unbound). Methods are normally evaluated by the percentage of hits, or approximately correct pre- dictions, commonly considered be those where the ligand RMSD error below 2Á [Alhossary et al., 2015; Hassan et al., 2017; McNutt et al., 2021].
Search-based docking methods. Traditional docking methods [TTout & Olson, 2010; Halgren al., 2004: Thomsen & Christensen, 2006] consist of a parameterized physics- scoring function and search algorithm. The scoring- funttiion takes in 3D structures and returns an estimate of the qual- ityllikelihood of the given pose, while the search stochastically modifies the ligand pose (position, orientation, and torsion angles) with the goal of finding the scoring function's 's global optimum. Re- cently, machine learning has been applied to parameterize the mooring-ffnction [McNutt et al., 2021; Méndez-Lucio Lucio et al., 2021]. J. However, these search-based methods remain computationally expen- sive to run, are often inaccurate when faced with the vast search space characterizing blind docking [Stärk et al., 2022], and significantly suffer when presented with apo-structures [Wong et al., 2022].
Machine learning for blind docking, Recently, EquiBind [Stärk et al., 2022] has tried to tackle the blind docking task by directly predicting pocket keypoints on both ligand and protein and aligning them. TANKBind [Lu al., 2022] improved over this by independently predicting docking pose (in the form of an interatomic distance matrix) for each possible pocket and then ranking them. Although these one-shot or few-shot regression-based prediction methods are orders of magnitude faster, their performance has not yet reached that of traditional search-based methods.
Diffusion generative models. Let the data distribution be the initial distribution px(x) of a con- tinuous diffusion process described by dxx = f(x,f) dt + g() dw, where w is the Wiener process. Diffusion generative models (DGMs) model the score² Vx log (x) of the diffusing data distribu- tion in order to generate data via the reverse diffusion dx [f(x f) g(t)² log Pt (x)] gt) dw [Song et al., 2021]. In this work, we always take f(x, t) 0. Several DGMs have been developed for molecular ML tasks, including molecule generation [Hoogeboom et al., 2022], conformer gen- eration [Xu et al., 2021], and protein design [Trippe et al., 2022]. However, these approaches learn distributions over the full Euclidean space R-3n with 3 coordinates per atom, making them ill-suited for molecular docking where the degrees of freedom are much more restricted (see Appendix E.2).
3 DOCKING AS GENERATIVE MODELING
Although Equi and other ML methods have provided strong runtime improvements by avoiding an expensive search process, their performance has not reached that of search-based methods. As our analysis below argues, this may be caused by the models' uncertainty and the optimization of an objective that does not correspond to how molecular docking is used and evaluated in practice.
Molecular docking objective. Molecular docking plays a critical role in drug discovery because the prediction of the 3D structure of bound ppoccin-ligane complex enables further computational and human expert analyses on the strength and properties of the binding interaction. Therefore, a docked prediction is only useful if its deviation from the true structure does not significantly affect the output of such analyses. Concretely, a prediction is considered acceptable when the distance between the structures (measured in terms of ligand RMSD) is below some small tolerance on the order of the length scale of atomic interactions (a few Ångström). Consequently, the standard evaluation metric used in the field has been the percentage of predictions with ligand RMSD (to the crystal ligand pose) below some value e.
Also known as diffusion probabilistic models (DPMs), denoising diffusion probabilistic models (DDPMs), diffusion models, score-hased models (SBMs), or score-based generative models (SGMs). Not to be confused with the scoring function of traditional docking methods.
3

---

Published as a conference paper at ICLR 2023
Crystal
EquiBind
TANKBind
DiffDock samples
DiffDock top-1
Figure 2: "DifFDocK top-1" refers the sample with the highest confidence. "DIFFDoCK sam- ples" to the other diffusion model samples. Left: Visual diagram of the advantage of generative models over regression models. Given uncertainty in the correct pose (represented by the orange distribution), regression models tend to predict the mean the distribution, which may lie a region of low density. Center: when there is a global symmetry in the protein (aleatoric uncer- tainty), EquiBind places the molecule in the center while DIFTDOCK is able to sample all the true poses. Right: even the absence of strong aleatoric uncertainty, the epistemic uncertainty causes EquiBind's prediction to have steric clashes and TANKBind's ' to have many self- eeff-intersections
However, the objective of maximizing the proportion of predictions with RMSD within some toler- ance e is not differentiable and cannot be used for training with stochastic gradient descent. Instead, maximizing the expected proportion of predictions with RMSD e corresponds to maximizing the likelihood of the true structure under the model's s output distribution, in the limit as goes to O. This observation motivates training generative model to minimize an upper bound on the negative log-likelihood of the observed structures under the model's ' distribution. Thus, we view molecu- lar docking as the problem of learning a distribution over ligand poses conditioned on the protein structure and develop diffusion generative model over this space (Section 4).
Confidence model. With a trained diffusion model, possible to sample an arbitrary number of ligand poses from the posterior distribution according to the model. However, researchers are often interested in seeing only one or a small number predicted poses and an associated confidence measure for downstream analysis. Thus, train confidence model the diffusion model and rank them based on its confidence that they are within the error tolerance. The top-ranked ligand pose and the associated confidence are then taken DIFFDocK's top-1 prediction and confidence score.
the poses sampled by
Problem with regression-basse methods. The difficulty with the development of deep learning models for molecular docking lies in the aleatoric (which is the data inherent uncertainty, e.g., the ligand might bind with multiple poses to the protein) and epistemic uncertainty (which arises from the complexity of the task compared with the limited model capacity and data available) on the pose. Therefore, given the available vo-variate information (only protein structure and ligand identity), any method will exhibit uncertainty about the correct binding pose among many viable alternatives. Any regression-style method that forced to select a single configuration that minimizes the expected square error would learn to predict the (weighted) mean of such alternatives. In contrast, generative model with the same co- varvaiit information would instead aim to capture the distribution over the alternatives, populating all/most of the significant modes even if similarly unable to distinguish the correct target. This behavior, illustrated in Figure 2, causes the regression- based models to produce significantly more physically implausible poses than method. In particular, we observe frequent steric clashes (e.g., 26% of EquiBind's predictions) preddcctiooooo and refffintersections in EquiBind's 's and TANKBind's predictions (Figures and 12). We found no intersections DIFFDOCK predictions. Visualizations and quantitative evidence of these phenomena are in Appendix F.I.
For example, the DDDT confidence score of AlphaFold2 [Jumper et al., 2021] has had a very significant impact in many applications [Necci al., 2021; Bennett al., 2022].
4

---

Published as a conference paper at ICLR 2023
4 METHOD
4.1 OVERVIEW
A ligand pose is an assignment of atomic positions in RR, so in principle, we can regard a pose x as an element in R where nt is the number of atoms. However, this encompasses far more degrees of freedom than are relevant in molecular docking. In particular, bond lengths, angles, and small rings in the ligand are essentially rigid, such that the ligand flexibility lies almost entirely in the torsion angles at rotatable bonds (see Appendix E.2 for further discussion). Traditional docking methods, as well as most ML ones, take as input seed conformation E R-3m of the ligand in isolation and change only the relative position and the torsion degrees of freedom in the final bound conformation." 4 The space of ligand poses consistent with is, therefore, an (m + 6) dimmennniaa submanifold Me C. L where m is the number of rotatable bonds, and the six additional degrees of freedom come from rototranslations relative to the fixed protein. We follow this paradigm of taking as input a seed conformation c, and formulate molecular docking as learning a probability distribution Pe (x y) over the manifold Me, conditioned on protein structure y.
DGMs on submanifolds have been formulated by De Bortoli et al. [2022] in terms of projecting diffusion in ambient space onto the submanifold. However, the kernel p(x; Xn) of such diffusion is not available in closed form and must be sampled numerically with a geodesic random walk, making training very inefficient. We instead define a one-to-one mapping to another, "nicer" manifold where the diffusion kernel can be sampled directly and develop a DGM in that manifold. To start, we restate the discussion the last paragraph as follows:
Any ligand pose consistent with seed conformation can be reached by combination of (1) ligand translations, (2) ligand rotations, and (3) changes to torsion angles.
This can be viewed as an informal definition of the manifold MM. Simultaneously, it suggests that given a continuous family of ligand pose transformations corresponding to the m + 6 degrees of freedom, distribution on Me can be lifted to a distribution on the product space of the correspond- ing groups
groups-which which is itself manifold. We will then show how to sample the diffusion kernel on this product space and train DGM over it.
4.2 LIGAND POSE TRANSFORMATIONS
We associate translations of ligand position with the 3D translation group T(3), rigid rotations of the ligand with the 3D rotation group S3(3), and changes in torsion angles at each rotatable bond with copy of the 2D rotation group (2)). More formally, we define operations of each of these groups on a ligand pose c E TR 3m The translation A : T(3) x 3 - RR is defined straightforwardly as Ag (r, x); = x; +r using the isomorphism T(3) = R3 where x; E R is the position of the ith atom. Similarly, the rotation A SO(3) x 3 + 3 is defined by Aro [R, x); = R(x: -x) + X where n Σ x, corresponding to rotations around the (unweighted) center of mass of the ligand.
Many valid definitions of a change in torsion angles are possible, as the torsion angle around any bond (a,, by) can be updated by rotating the a side, the 6, side, or both. However, we can specify changes of torsion angles to be disentangled from rotations or translations. To this end, we define the operation of elements of 50()) such that it causes a minimal perturbation (in an RMSD sense) to the structure:?
Let BN,0. (x) E 3m be any valid torsion update by 02 around the kth rotatable bond (ax, bx). We define A (2( x 3g + 3 such that
Definition.
A 0 x) RMSDAlignxx, (B1,0, o... Bm, mm. (x))
where 0 = (01,... .m)) and RMSDAlign(X,xx) x (xx
arg min
RMSD(x, x!)
(1)
'RDKit ETKDG is popular method for predicting the seed conformation. Although the structures may not be predicted perfectly, the emrrr largely the torsion angles, which are resampled anyways. SSince we do not define or use the composition of elements of SO(())" strictly speaking, it is a product space but not group and can be alternatively thought of as the torus TTT with an origin element.
5

---

Published as a conference paper at ICLR 2023
This means that apply all the m torsion updates in any order and then perform global RMSD alignment with the unmodified pose. The definition is motivated by ensuring that the infinitesimal effect of torsion is orthogonal to any rototranslation, i..., it induces no linear angular momen- tum. These properties can be stated more formally as follows (proof in Appendix A):
Proposition I. Let y(t) A x) for some O and where 10 = (101,.... -t0m). Then the linear and angular momentum are zero: at
d V (t= 0 and (x x) Viloo 0 where x 1 E
Now consider the product space" P = T x 50(3) x S0( TM and define A IP x 3
R as
A((r, R 0), x) Ar(r, A (R, A (0, x)))
(2)
These definitions collectively provide the sought-after product space corresponding to the docking degrees of freedom. Indeed, for seed ligand conformation c, we can formally define the space of ligand poses M. {g c) E P.. This corresponds precisely to the intuitive notion of the space of ligand poses that can be reached by rigid- body motion plus torsion angle flexibility.
4.3 DIFFUSION ON THE PRODUCT SPACE
We now proceed to show how the product space can be used to learn a DGM over ligand poses in Me. First, we need theoretical result (proof in Appendix A): Proposition 2. For a given seed conformation c, the map A(t, : Me is a bijection.
which means that the inverse A-1 : Me - IP given by A(g, c) + 9 maps ligand poses x E Me to points on the product space IP. We are now ready to develop a diffusion process on P.
De Bortoli et al. [2022] established that the DGM framework transfers straightforwardly to Rie- mannian manifolds with the score and score model as elements of the tangent space and with the geodesic random walk the reverse SDE solver. Further, the score model can be trained in the standard manner with denoising score matching [Song & Ermon, 2019]. Thus, to implement a dif- fusion model on P., it suffices develop method for sampling from and computing the score of the diffusion kernel on P. Furthermore, since IP is a product manifold, the forward diffusion pro- ceeds independently in each manifold [Rodolà et al., 2019], and the tangent space is direct sum: T.P = T3 @ TR (3) @ (2) = R + R + R where = (r, R, 0). Thus, it suffices to sample from the diffusion kernel and regress against its score in each group independently.
In all three groups, we define the forward SDE dx V (t)/ dw where o² oie of of for T(3), S3(3), and SO(2) respectively and where is the corresponding Brownian motion. Since T(3) = IRa the translational case is trivial and involves sampling and computing the score of standard Gaussian with variance o2 ot). The diffusion kernel on SO(3) is given by the IGSO(() distribution [Nikolayev & Savyolov, 1970; Leach et al., 2022], which can be sampled in the axis-angle parameterization by sampling unit vector @ E so(3) uniformly? and random angle E [0, π according to
p(w) - cos f(w) where f(w) (00=>(e (21 + 1) exxe(p((((( 1)o2 /2) sin + 1/22)w)
00
(3)
sin(w/2)
π
Further, the score of the diffusion kernel Il R) = (4 log f(w)) E TR SO(3),
where R' = R(wx)) R is the result of applying Euler vector we to R YYmm et al., 2023]. The score computation and sampling can be accomplished efficiently by precomputing the truncated infinite series and interpolating the CDF of p(w), respectively. Finally, the SO(2)T" group is diffeomorphic to the torus T" on which the diffusion kernel is wrapped normal distribution with variance (t).. This can be sampled directly, and the score can be precomputed as a truncated infinite series [Jing et al., 2022].
4.4 TRAINING AND INFERENCE
Diffusion model. Although have defined the diffusion kernel and score matching objectives on P. we nevertheless develop the training and inference procedures to operate on ligand poses in 3D
6 Since we never compose elements of P, we do not need to define a group structure. so()) is the tangent space of SO(3) at the identity and is the space of Fuler (or rotation) vectors, which equivalent the axis-angle parameterization.

---

Published as a conference paper at ICLR 2023
coordinates directly. Providing the full 3D structure, rather than abstract elements of the product space, to the score model allows it to reason about physical interactions using SE()) equivariant models, not be dependent on arbitrary definitions of torsion angles [Jing al., 2022], and better generalize to unseen complexes. In Appendix B, we present the training and inference procedures and discuss how resolve their dependence on the choice of seed conformation used to define the mapping between Me and the product space.
Confidence model. In order to collect training data for the confidence model d(x, y), we run the trained diffusion model to obtain a set of candidate poses for every training example and generate labels by testing whether or not each pose has RMSD below 2Å. The confidence model then trained with cross-entropy loss to correctly predict the binary label for each pose. During inference, the diffusion model is run to generate N poses in parallel, which are passed to the confidence model that ranks them based on its confidence that they have RMSD below 2Å.
4.5 MODEL ARCHITECTURE
We construct the score model s(x,y, 1) and the confidence model d(x,y) to take as input the current ligand pose x and protein structure y in 3D space. The output of the confidence model is a single scalar that is SE(3)-invariant (with respect to joint rototranslations of x,y) as lig- and pose distributions are defined relative to the protein structure, which can have arbitrary lo- cation and orientation. On the other hand, the output of the score model must be in the tangent space T @ (3)) @ (o)O(2) The space T- 2 T corresponds to translation vectors and (3)) = R3 to rotation (Euler) vectors, both of which are -e)eeqqquivariint. Finally, (o2(()) corresponds to scores on (3))iinnrriia quantities (torsion angles). Thus, the score model must predict two (EE())-quivarri vectors for the ligand whole and an S (iin)-iinari scalar at each of the m freely rotatable bonds.
The score model and confidence model have similar architectures based on (qqiivravvv con- volutional networks over point clouds [Thomas et al., 2018; Geiger et al., 2020]. However, the score model operates the confidence model has access to the all-atom structure. This multiscale setup yields improved performance and significant speed-up w.r.t. doing the whole process at the atomic scale. The architectural components are summarized below and detailed in Appendix C.
coarse-grained representation of the protein with a-carbon atoms, while
Structures are represented as heterogeneous geometric graphs formed by ligand atoms, protein residues, and (for the confidence model) protein atoms. Residue nodes receive as initial features language model embeddings trained on protein sequences [Lin et al., 2022]. Nodes are sparsely connected based on distance cutoffs that depend on the types of nodes being linked and on the diffu- sion time. The ligand atom representations after the final interaction layer are then used to produce the different outputs. To produce the two IR3 vectors representing the translational and rotational scores, we convolve the node representations with a tensor product filter placed at the center of mass. For the torsional score, we use pseudotorque convolution to obtain scalar at each rotatable bond of the ligand analogously to Jing et al. [2022], with the distinction that, since the score model operates coarse -griined representations, the output is not a pseudoscalar (its parity is neither odd nor even). For the confidence model, the single scalar output is produced by mean- meean-poling the ligand atoms' scalar representations followed by fully connected layer.
5 EXPERIMENTS
Experimental setup. We evaluate our method on the complexes from PDBBind [Liu et al., 2017], a large collection of protein-ligand structures collected from PDB [Berman et al., 2003], which was used with time-based splits to benchmark many previous works [Stärk et al., 2022: Volkov et al., 2022; Lu et al., 2022]. We compare DIFFDOCK with state-of-the- art search-based methods SMINA [Koes et al., 2013]. Quick [Hassan et al., 2017], GLIDE [Halgren et al., 2004], and GNINA [McNutt et al., 2021] as well as the older Autodock Vina [Trott& & Olson, 2010], and the recent deep learning methods EquiBind and TANKBind presented above. Extensive details about the experimental setup, data, baselines, and implementation available at ht tps : /ttp:p::///ttitte //ghtubb. con/ htt///coroooooooooo /D/i/fif DDokk. The repository also contains videos of the reverse diffusion process (images of the same are in Figure 14).
in Appendix D.4 and all code is
7

---

Published as a conference paper at ICLR 2023
Table 1: PDBBind blind docking, All methods receive small molecule and are tasked find its binding location, orientation, and conformation. Shown is the percentage of predictions with RMSD < 2Á and the median RMSD (see Appendix D.3). The top half contains methods that directly find the pose; the bottom half those that use pocket prediction method. The last two lines show our method's performance. In parenthesis, we specify the number of poses sampled from the generative model. indicates that the method runs exclusively on CPU, "-" means not applicable; some cells are empty due to infrastructure constraints. For TANKBind, the runtimes for the top-1 and top-5 predictions are different. No method received further training or tuming ESMFold structures. Further evaluation details are in Appendix D.4 and further metrics are reported in Appendix F.
Holo crystal proteins
Apo ESMFold proteins
Top-1
RMSD
Top-5
RMSD
Top-1
RMSD
Top-5
RMSD
Average
Method
%<2
Med.
%<2
Med.
%<2
Med.
%<2
Med.
Runtime (s)
GNINA
22.9
7.7
32.9
4.5
2.0
22.3
4.0
14.22
127
SMINA
18.7
7.1
29.3
4,6
3.4
15.4
6.9
10.0
126
GLIDE
21.8
9.3
1405
EQUIBIND
5.5
6.2
1.7
7.1
0.04
TANKBIND
20.4
4.0
24.5
3.4
10.4
5.4
14.7
4.3
0.7/2.5
SMRANK+SMINA
20.4
33.2
4.4
4.6
10.0
10.3
7.0
126*
P2RANK+GNINA
28.8
5.5
38.3
3.4
8.6
11.2
12.8
7.2
127
EQUIBIND+SMINA
23,2
6.5
38.6
3,4
4.3
8.3
11.7
5.8
126
EQUIBIND+GNINA
28.8
4.9
39.1
3.1
10.2
8.8
18.6
5.6
127
DIFFDOCK (10)
35.0
3.6
40.7
2.65
21.7
5.0
31.9
3.3
10
DIFFDOCK (40)
38.2
3.3
44.7
2.40
20.3
5.1
31.3
3.3
40
As we evaluating blind docking, the methods receive two inputs: the ligand with predicted seed conformation (e.g., from RDKit) and the crystal structure of the protein. Since search-based methods work best when given starting binding pocket restrict the search space, we also test the combination of using an ML-based method, such as P2Rank [Krivák & Hoksza, 2018] (also used by TANKBind) or EquiBind to find an initial binding pocket, followed by search-based method to predict the exact pose in the pocket. To evaluate the generated poses, we compute the heavy-atom RMSD (permutation symmetry corrected) between the predicted and the ground-truth ligand when the protein structures are aligned. Since all methods except for EquiBind generate multiple ranked structures, we report the metrics for the highest-ranked prediction as the top-1 prediction and top-5 refers to selecting the most accurate pose out of the highest ranked predictions; useful metric when multiple predictions are used for downstream tasks.
Apo-structure docking setup. In order to test the performance of the models against computa- tionally generated apo-structures, we run ESMFold [Lin et al., 2022] on the proteins of the test-set of PDBBind and align its predictions with the crystal holo-structure. In order to align structure changing as little as possible to the conformation of the binding site, we use in the alignment of the residues a weight exponential with the distance to the ligand. More details on this alignment are provided in Appendix D.2. The docking methods are then run against the generated structures and compared with the ligand pose inferred by the alignment.
Docking accuracy. DIFFDOCK significantly outperforms all previous methods (Table 1). In par- ticular, DIFFDOCK obtains an impressive 38.2% top-1 success rate (i.e.. percentage of predictions with RMSD <2 when sampling 40 poses and 35.0% when sampling just 10. This performance vastly surpasses that of -fate-of-the-are commercial software such as GLIDE (21.8%. 2=2 10- and the previous trat-oofthh--a deep learning method TANKBind (20.4%, of ML-based pocket prediction in combination with search-based docking methods improves over the baseline performances, but even the best of these (EquiBind+GNINAA reaches a success rate of only 88.8% (p=0.0003). Unlike regression methods like EquiBind, DIFFDOCK is able to provide multiple diverse predictions of different likely poses, highlighted in the top-5 performances.
The use
Apo-structure docking accuracy. Previous work Wong et al. [2022], that highlighted that existing docking methods are not suited to dock to (computational) apo-structures, validated by the results
Most commonly used evaluation metric [Alhossary et al., 2015; Hassan al., 2017; McNutt al., 2021]
8

---

Published as a conference paper at ICLR 2023
Fraction with RMSD < <2
Percentage w/ RMSD -2
0.45
0.35
Top-1 performance
40
Top performance
Baseline performance
0.25
Top-10 performance
20
Perfect selection
Perfect selection
Confidence model
30
60
Number generative samples
Percentage of rejected complexes
Figure 3: Left: DIFFDOCK's performance as a function of the number of samples from the gen- erative model. "Perfect selection" refers choosing the sample with the lowest RMSD. Right: Selective accuracy: percentage of predictions with RMSD below 2Å when only making predictions for the portion of the dataset where DIFFDOCK is confident.
in Table where previous methods obtain top-1 accuracies of only 10% or below. This is most likely due to their reliance on trying to find key-lock matches, which makes them inflexible to imperfect protein structures. Even when built-in options allowing side-chain flexibility are activated, the results do not improve (see Appendix 8). Meanwhile, DIFFDOCK is able to retain larger proportion of its accuracy, placing the top-ranked ligand below 2 away on 22% of the complexes. imperfect structures, even without retraining, can be attributed combination of (1) the robustness of the diffusion model to small perturbations in the backbone atoms, and (2) the fact that DIFFDOCK does not use the exact position of side chains the score model and is therefore forced to implicitly model their flexibility.
This ability to generalize
Inference runtime. DIFFDOCK holds its superior accuracy while being (on GPU) to 12 times faster than the best search-based method, GNINA (Table 1). This high speed is critical for applica- tions such as high throughput virtual screening for drug candidates or reverse screening for protein targets, where one often searches over vast number of complexes. As diffusion model, DIFF- DOCK is inevitably slower than the one-shot deep learning method EQUIBIND, Figure 3-left and Appendix F.3, it can be significantly sped up without significant loss of accuracy.
as shown in
Selective accuracy of confidence score. As the top-1 results show, DIFFDOCK'S ' confidence model very accurate in ranking the sampled poses for given complex and picking the best one. We also investigate the selective accuracy of the confidence model across different complexes by evaluating how DIFFDOCK's 's accuracy increases if it only makes predictions when the confidence is above certain threshold, known as selective prediction. In Figure 3-right, we plot the success rate as we decrease the percentage of complexes for which we make predictions, i.e., increase the confidence threshold. When only making predictions for the top one-third of complexes in terms of model confidence, the success rate improves from 38% 83%. Additionally, there a Spearman correla- tion of 0.68 between DIFFDoCK's 's confidence and the negative RMSD. Thus, the confidence score good indicator of the quality of DIFFDOcK's top-ranked sampled pose and provides highly valuable confidence measure for downstream applications.
CONCLUSION
We presented DIFFDOCK, diffusion generative model tailored the task of molecular docking. This represents a paradigm shift from previous deep learning approaches, which use regression- based frameworks, molecular docking. To produce a fast and accurate generative model, we designed a diffusion pro- cess over the manifold describing the main degrees of freedom of the task via ligand pose transfor- mations spanning the manifold. Empirically, DIFFDOCK outperforms the state-of-the-art by large margins on PDBBind, has fast inference times, and provides confidence estimates with high selective accuracy. Moreover, unlike previous methods, it retains large proportion of its accuracy when run on computationally folded protein structures. Thus, DIFFDOCK can offer great value for many exist- ing real-world pipelines and opens up new avenues of research on how to best integrate downstream tasks, such as affinity prediction, into the framework and apply similar ideas to protein-protedd and protein-nucleie acid docking.
generative modeling approach that is better aligned with the objective of

---

Published as a conference paper at ICLR 2023
ACKNOWLEDGMENTS
We pay tribute to Octavian-Bugen Ganea (1987-2022), dear colleague, mentor, and friend without whom this work would have never been possible.
We thank Rachel Wu, Jeremy Wohlwend, Felix Faltings, Jason Yim, Victor Quach, Saro Passaro, Patrick Walters, Michael Heinzinger, Mario Geiger, Michael John Arcidiacono, Noah Getz, and John Yang for valuable feedback and insightful discussions. We thank Wei Lu for his help with running TANKBind. We thank Valentin De Bortoli, Emile Mathieu, Brian Trippe, and Jason Yim for the critical discussions and help formalize diffusion 50(3). Noah Getz contributed the implementation of the apo-structures alignment procedure.
This work supported by the Machine Learning for Pharmaceutical Discovery and Synthesis (MLPDS) consortium, the Abdul Latif Jameel Clinic for Machine Learning in Health, the DTRA Discovery of Medical Countermeasures Against New and Emerging (DOMANE) threats program, the DARPA Accelcrated Molecular Discovery program and the Sanofi Computational Antibody Design grant. Bowen Jing acknowledges the support from the Department of Energy Computational Science Graduate Fellowship.
REFERENCES
Amr Alhossary, Stephanus Daniel Handoko, Yuguang Mu, and Chee-Keong Kwoh, Fast, accurate, and reliable molecular docking with QuickVina 2. Bioinformatics, 31(13):2214-23 2216. 02 2015.
Nathaniel Bennett, Brian Coventry, Inna Goreshnik, Buwei Huang, Aza Allen, Dionne Vafeados, Ying Po Peng, Justas Dauparas, Minkyung Back, Lance Stewart, et al. Improving de novo protein binder design with deep learning. bioRxiv, 2022.
H. Berman, K. Henrick, and H. Nakamura. Announcing the worldwide Protein Data Bank. Nat Struct Biol, 10(12)9980, Dec 2003.
Richard H Byrd, Peihuang Lu, Jorge Nocedal, and Ciyou Zhu. A limited memory algorithm for bound constrained optimization. SIAM Journal on scientific computing, 1995.
Valentin De Bortoli, Emile Mathieu, Michael Hutchinson, James Thornton, Yee Whye Teh, and Arnaud Doucet. Riemannian score-based generative modeling. arXiv preprint 32X76, 2022.
Mario Geiger and Tess Smidt. c3nn: Euclidean neural networks. arXiv preprint arXi::2207.09453, 2022.
Mario Geiger, Tess Smidt, Alby M., Benjamin Kurt Miller, Wouter Boomsma, Bradley Dice, Kos- tiantyn Lapchevskyi, Maurice Weiler, Michat Tyszkiewicz, Simon Batzner, Martin Uhrin, Jes Frellsen, Nuri Jung, Sophia Sanbor,, Josh Rackers, and Michael Bailey. Euclidean neural net- works: e3mn, 2020.
Thomas A Halgren, Robert B Murphy, Richard A Friesner, Hege S Beard, Leah L Frye, W Thomas Pollard, and Jay L Banks. Glide: new approach for rapid, accurate docking and scoring. 2. enrichment factors in database screening. Journal of medicinal chemistry, 2004.
Nafisa M. Hassan, Amr A. Alhossary, Yuguang Mu, and Chee- Keong Kwoh. Protein-ligand blind docking using quickvina-w with inter-process spatio-temporal integration. Scientific Reports, (1):15451, Nov 2017.
Emiel Hoogeboom, Victor Garcia Satorras, Clement Vignac, and Max Welling. Equivariant dif- fusion for molecule generation in 3d. In International Conference Machine Learning, pp. 8867-8887. PMLR, 2022.
Bowen Jing, Gabriele Corso, Jeffrey Chang, Regina Barzilay, and Tommi Jaakkola. Torsional diffu- sion for molecular conformer generation. arXiv preprint arXiv:2206.01729, 2022.
10

---

Published as a conference paper at ICLR 2023
John Jumper, Richard Evans, Alexander Pritzel, Tim Green, Michael Figurnov, Olaf Ronneberger, Kathryn Tunyasuvunakool, Russ Bates, Augustin Zidek, Anna Potapenko, et al. Highly accurate protein structure prediction with alphafold. Nature, 596(7873):583-58) 2021.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv: (ee..1412.660 2014.
David Ryan Koes, Matthew P Baumgartner, and Carlos J Camacho. Lessons learned in empirical scoring with smina from the esar 2011 benchmarking exercise. Journal of chemical information and modeling, 13(8):1893-1900) 1993-1904, 2013.
Radoslav Krivák and David Hoksza. P2rank: machine learning based tool for rapid and accurate prediction of ligand binding sites from protein structure. Journal of cheminformatics, 2018.
Adam Leach, Sebastian M Schmon, Matteo T Degiacomi, and Chris G Willcocks. Denoising diffu- sion probabilistic models on so (3) for rotational alignment. In ICLR 2022 Workshop on Geomet- rical and Topological Representation Learning, 2022.
Zeming Lin, Halil Akin, Roshan Rao, Brian Hie, Zhongkai Zhu, Wenting Lu, Allan dos San- tos Costa, Maryam Fazel-Zarandi, Tom Sercu, Sal Candido, and Alexander Rives. Language models of protein sequences at the scale evolution enable accurate structure prediction. arXiv, 2022.
Zhihai Liu, Minyi Su, Li Han, Jie Liu, Qifan Yang. Yan Li, and Renxiao Wang. Forging the basis for developing pootiin-lggane interaction scoring functions. Accounts of Chemical Research, 50 (2):302-309 2017.
Wei Lu, Qifeng Wu, Jixian Zhang, Jiahua Rao, Chengtao Li, and Shuangjia Zheng. Tankbind: Trigonometry neural networks for drug-protein binding structure prediction. Advances in neural information processing systems, 2022.
Andrew T McNutt, Paul Francoeur, Rishal Aggarwal, Tomohide Masuda, Rocco Meli, Matthew Ragoza, Jocelyn Sunseri, and David Ryan Koes. Gnina 1.0: molecular docking with deep learn ing. Journal of cheminformatics, 13(1):1-20,202
Rocco Meli and Philip C. Biggin. spyrmsd: cymmetry-corrected rmsd calculations in python. Jour- nal of Cheminformatics, 12(1):49, 2020.
Oscar Méndez-Lucio, Mazen Ahmad, Ehecatl Antonio del Rio- Rio-Chanona, and Jörg Kurt Wegner. A geometric deep learning approach to predict binding conformations of bioactive molecules. Nature Machine Intelligence, 10333)003-1 2021.
Marco Necci, Damiano Piovesan, and Silvio CE Tosatto. Critical assessment of protein intrinsic disorder prediction. Nature methods, 18(5):472-481, 2021.
Dmitry Nikolayev and Tatjana I Savyolov. Normal distribution on the rotation group so (3). Tex- fures and Microstructures, 29, 1970.
Srinivas Ramachandran, Pradeep Kota, Feng Ding, and Nikolay V Dokholyan. Automated mini- mization of steric clashes in protein structures. Proteins, 79(1):261-277 January 2011.
Emanuele Rodolà, Zorah Lähner, Alexander M Bronstein, Michael M Bronstein, and Justin Solomon. Functional maps representation product manifolds. In Computer Graphics Forum, volume 38, pp. 678-689. Wiley Online Library, 2019.
Kristof Schütt, Pieter-Jan Kindermans, Huziel Enoc Sauceda Felix, Stefan Chmiela, Alexandre Tkatchenko, and Klaus-Robert Müller. Schnet: continuous-filte convolutional neural network for modeling quantum interactions. Advances in neural information processing systems, 2017.
Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. Advances in Neural Information Processing Systems, 32, 2019.
11

---

Published as a conference paper at ICLR 2023
Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. In Interna- tional Conference on Learning Representations, 2021.
Hannes Stärk, Octavian Ganea, Lagnajit Pattanaik, Regina Barzilay, and Tommi Jaakkola. Equibind: Geometric deep learning for drug binding structure prediction. In International Conference on Machine Learning. pp. 20503-205211 PMLR. 2022.
Nathaniel Thomas, Tess Smidt, Steven Kearnes, Lusann Yang, Li Li, Kai Kohlhoff, and Patrick Riley. Tensor field networks: Rotation-and translation-equiiariant neural networks for 3d point clouds. arXiv preprint, 2018.
René Thomsen and Mikael H Christensen. Moldock: a new technique for high-accuracy molecular docking. Journal of medicinal chemistry, 49(1 (9(11):331-3321, 2006.
Brian Trippe, Jason Yim, Doug Tischer, Tamara Broderick, David Baker, Regina Barzilay, and Tommi Jaakkola. Diffusion probabilistic modeling of protein backbones in 3d for the motif- scaffolding problem. arXiv preprint arXiv:2206. 04. 19, 2022.
Oleg Trott and Arthur Olson. Autodock vina: improving the speed and accuracy of docking with new scoring function, efficient optimization, and multithreading. Journal of computational chemistry, 31(2)455-461, 2010.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, tukasz Kaiser, and Illia Polosukhin. Attention is all you need. Advances in neural informa- tion processing systems, 2017.
Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Courna- peau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nel- son, Eric Jones, Robert Kern, Eric Larson, Carey, Îlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mul- bregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 2020a.
Pauli Virtanen, Ralf Gommers, Travis E Oliphant, Matt Haberland, Tyler Reddy, David Courna- peau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, al. Scipy 1.0: fundamental algorithms for scientific computing in python. Nature methods, 2020b.
Mikhail Volkov, Joseph-André Turk, Nicolas Drizard, Nicolas Martin, Brice Hoffmann, Yann Gaston-Mathé, and Didier Rognan. On the frustration to predict binding affinities from protein- ligand structures with deep neural networks. Journal of Medicinal Chemistry, 65(11);7946-7958 Jun 2022. ISSN 0022-2623.
Felix Wong, Aarti Krishnan, Erica Zheng, Hannes Stärk, Abigail I. Manson, Ashlee Earl, Tommi Jaakkola, and James Collins. Benchmarking alphafold-enabled molecular docking pre- dictions for antibiotic discovery. Molecular Systems Biology, 2022.
Minkai Xu, Lantao Yu, Yang Song, Chence Shi, Stefano Ermon, and Jian Tang. Geodiff: A geo- metric diffusion model for molecular conformation generation. In International Conference on Learning Representations, 2021.
Jason Yim, Brian L Trippe. Valentin De Bortoli, Emile Mathieu, Arnad Doucet, Regina Barzilay, and Tommi Jaakkola. Se(3) diffusion model with application to protein backbone generation. arXiv preprint, 2023.
12

---

Published as a conference paper at ICLR 2023
A PROOFS
A.1 PROOF OF PROPOSITION
Proposition 1. Let yt) (t0, x) for some I and where to (#01,... t0mm. Then the linear and angular momentum are zero: d (1=0 and ( x) 2: Vit where x 1/2 -i X
Proof. Let y(t) = (Bt) 0, x) - x) + p(t) where Bt,, 0 = Bi, .... ....Bm.t.... and R(t), p()) the rotation (around x) and translation associated with the optimal RMSD alignment between Bt, 0, x) and x. By definition of RMSD, for any 1, R(t) and p(t) minimize
I|y(t) -x (B(( 0, x) + pt) xxl
(4)
For infinitesimal dt, the RHS becomes
RHS = =[ 0, x) + p(dl) - x||
= =|[ [R (0) dt + R())) (B, dt + (00, 0, - x) + x + p'(0) dt + pp() - x
(5)
= (0) (x - x) + 0,) +p((0) || dl
where we have used R[(0) [(0, 0,x) xx. and p(0) 0. Thus, see that RMSD alignment implies that the derivatives of P(()).(() minimize the norm of
y'0)) R' (x)
(0 +p((0)
This expression represents the instantancous velocity of the points y, at t = 0. We now show that minimizing the velocity results in zero linear and angular momentum.
We abbreviate (t,, 0, x())); == b, and p' = v. Further, let r; = x; - &, such that the rotational contribution to the velocity can be written in terms of angular velocity vector With this, at t t= we have
y/ +
We thus obtain the squared norm as
E = - +wxr++v (
Ellb:ll 2b. (wxr)) 2b. (w x..) + [[v]]]
-bb 2w E(rxb:)) +2 Σbb, +
(8)
where we have used the fact that Σ., 0 and where I(r) x inertia tensor. To minimize the squared norm (and thus the norm itself), we set gradients with respect to v, w to zero. This gives
-. is the
Σ b. and w w- -I()--
Σ ri x
(9)
Now with y/ b, r we evaluate the linear momentum
1 Σxー= Σb. + Eri+ + =
(10)
which is zero by direct substitution of v. Similarly, we evaluate the angular momentum
Σ r, x y = ΣEr r x + Σ r x (w (x)) + Σ r xv
(11)
Er,bb + I(r)
which is zero by direct substitution of w. Thus, the linear and angular momentum are zero at f = for arbitrary x.
Note that since did not use the particular form of B(t0, x) in the above proof, have shown that RMSD alignment can be used to disentangle rotations and translations from the infinitesimal action of any arbitrary function.
13

---

Published as a conference paper at ICLR 2023
A.2 PROOF OF PROPOSITION 2
Proposition 2. For a given seed conformation c, the map A(,, c) PP + Me is a bijection.
Proof. Since we defined Me = {A(g, c) g E P), A(,c) is automatically surjective. We now show that it is injective. Assume for the sake of contradiction that c(:,c) is not injective, so that there exist elements of the product space 91, 92 E P with 91 + 92 That is,
but with A(g1,c) = A(g2,
AI Ame Ami A ( c))) = (r, A R2, Aur (62.
(12)
which we abbreviate c(l) c(2) Since only Arr changes the center of mass E, ci/n, we have E, of /n = Σ, c,/n + 1ı and Σ, c,c/9//1 /n = Σ, c./n + r2. However, since c( = c(2) this implies 1ı = r2. Next, consider the torsion angles TI to some choice of dihedral angles at each rotatable bond. Because Ar and A are rigid-body motions, only Aur changes the dihedral angles; in particular, by definition we have 7.1) = 7i 0(1)
". of c(2) corresponding
mod 2T and T: (2) Tj + 0(2) mod 2% for all i = m. However, because T; 7
this means 8(1) = 8!2) for all and therefore 01 02 (as elements of SO(2))). Now denote
c = Aur (01, c) = A (02. c) and apply Au( ( = A-(-22.)) to both sides of Equation 12.
We then have
A R1, c²) Aroe R22, c*)
(13)
which further leads to
c' = R -(²
(14)
In general, this does not imply that R₁ Rg. However, R₁ R2 is possible only if is degenerate,
in the sense that all points are collinear along the shared axis of rotation of R1, Rg. However, in practice, conformers never consist of collinear set of points, so we can safely assume R₁ R2.
We now have (r1, R1, 01) = (r2, R2, 02), or 91 92, contradicting our initial assumption. We thus conclude that A(, c) is injective, completing the proof.
B TRAINING AND INFERENCE
In this section we present the training and inference procedures of the diffusion generative model. First, however, there are few subtleties of the generative approach to molecular docking that are worth mentioning. Unlike the standard generative modeling setting where the dataset consists of many samples drawn from the data distribution, each training example (x",y) of protein structure y and ground-truth ligand pose x" is the only sample from the corresponding conditional distribution Px y) defined over Mx* Thus, the innermost training loop iterates distinct conditional distributions Px*(( (- y), along with a single sample from that distribution, rather than over samples from a common data distribution Panta (x).
As discussed in Section 4, during inference, e is the ligand structure generated with a method such RDKit. However, during training we require Me Mx in order define a bijection between E <Mx' and P. If we take с E Mxx, there will be a distribution shift between the manifolds Me considered training time and those considered at inference time. To circumvent this issue, at training time we predict e with RDKit and replace x* with arg milxttM,, RMSD(x", x!) using the conformer matching procedure described Jing et al. [2022].
The above paragraph may be rephrased more intuitively as follows: during inference, the generative model docks a ligand structure generated by RDKit, keeping its non-torsional degrees of freedom (e.g., local structures) fixed. At training time, however, if we train the score model with the local structures of the ground truth pose, this will not correspond to the local structures seen at inference time. Thus, at training time, replace the ground truth pose by generating a ligand structure with RDKit and aligning to the ground truth pose while keeping the local structures fixed.
With these preliminaries, we now continue to the full procedures (Algorithms and 2). The training and inference procedures of a score-based diffusion generative model on a Riemannian manifold consist of (1) sampling and regressing against the score of the diffusion kernel during training; and (2) sampling a geodesic random walk with the score as drift term during inference [De Bortoli al., 2022]. Because we have developed the diffusion process on PP but continue to provide the score
14

---

Published as a conference paper at ICLR 2023
model with elements in Me C R the full training and inference procedures involve repeatedly interconverting between the two spaces using the bijection given by the seed conformation c.
Algorithm 1: Training procedure (single epoch)
Input:
Training pairs {(x*, RDKit predictions {c}
foreach c, x do
Let Xg +- arg min RMSD(xx,x)))
Compute (ro. 0g) < A- (xo);
Sample ~ Uni([O, 1]):
Sample Ar, A R2 A0 from diffusion kernels p. (), o) (-10) 10);
Set
+ Ar;
Set R2 + (AR)
Set O2 - 0o A0 mod 2m;
Compute x, + A((r, R4 0.) c):
Predict scores a E R23 3 E R3, Y E I = s(x,,c,y,
Take optimization step on loss
C a
log p Ar 0)1 ||
V logp
[)] 11 V Vlog p (A0 O)]
Algorithm 2: Inference procedure
Input: RDKit prediction c, protein structure y (both centered at origin)
Output: Sampled ligand pose X
Sample ON ~ Uni( (2) Rx Uni(SO(3)), rN ~ (O (T))
Let XN A(((N,, Ry 0)),
for π + N to 1 do
Let 1 = n/N and Do = oz (n/N) - oz((n - 1)/N) and similarly for Ao Ao 1Rm + s(x,, c, t);
tor
Predict scores a E 12,8 R3
Sample Zur, 2m,, Zior from NO, N,,NN ) Ao respectively;
Set 12-11 + In + Ao + Zy:
Set Rn-1 + - + Zr0e Ri:
Set 0n-1 1.+,,, 0n + (Ao tor 7 + Z mod 2m:
Compute Xn-1 A((rn-1, R 0 1n-1), c);
Return No:
However, as noted the main text, the dependence these procedures on the exact choice is potentially problematic, it suggests that at inference time, the model distribution may be different depending on the orientation and torsion angles of c. Simply removing the dependence of the score model on is not sufficient since the update steps themselves still occur P and require choice of c to be mapped to Me. However, notice that the update steps-in both training and inference- consist of (I) sampling the diffusion kernels at the origin; (2) applying these updates to the point on F; and (3) transferring the point on P to Me via c). Might it instead be possible to apply the updates to 3D ligand poses x E Me directly?
It turns out that the notion of applying these steps to ligand poses "directly" corresponds to the formal notion of group action. The operations Aq, A A that we have already defined are formally group actions if they satisfy (01 x) = A A(g2, x)). While true for Ar. Ana this is not generally true for A if we take ( to be the direct product group; however, the approximation is increasingly good the magnitude of the torsion angle updates decreases. If we then define IP be the direct product group of its constituent groups, A is group action of P on Me. as the operations of Ag. Acot A commute and are (under the approximation) individually group actions.
The implication of A being group action can be seen as follows. Let d 9a9. be the update which brings Jo E TP to gb E P via left multiplication, and let x,, X the corresponding ligand poses A(ga, c), A(g,, c). Then
X A(goga 9a. c) A(8,xa)
(15)
15

---

Published as a conference paper at ICLR 2023
which means that the updates ò can be applied directly to xa using the operation A. The training and inference procedures then become Algorithm and below. The initial conformer с is no longer used, except in the initial steps to define the manifold to find the closest point to x" in training, and to sample x from the prior over Me in inference.
Conceptually speaking, this procedure corresponds to "forgetting" the location of the origin element on Me. which is permissible because a change of the origin to some equivalent seed E MM merely translates-via right multiplication by - original and diffused data distributions on IP, but does not cause any changes on Me itself. The training and inference routines involve updates-formald formally left multiplications- corresponds to group actions on Mc. the updates can act on Me directly, without referencing the origin c.
group elements, but as left multiplication on the group
We find that the approximation of A group action works quite well in practice and use Algo- rithms and for all training and experiments discussed in the paper. Of course, disentangling the torsion updates from rotations in way that makes Aior exactly group action would justify the procedure further, and we regard this as a possible direction for future work.
Algorithm 3: Approximate training procedure (single epoch)
Input:
Training pairs {(x* y)), RDKit predictions {c)
foreach c, xx, do
Let Xg + arg min RMSD(xx",
Sample ~ Uni([O, 1]);
Sample Ar, AR, A0 from diffusion kernels y/(( (t, 0), 0);
Compute x, f A((Ar AR, x0):
Predict scores a E _², 3 E R3 y E IR s(x,, y,t)
Take optimization step loss
L
|o
V log p.
Ar 1] 2 ||3 v log 0)] 1b V Vlog (A0 0)]
Algorithm 4: Approximate inference procedure
Input: RDKit prediction c, protein structure y (both centered at origin)
Output: Sampled ligand pose X
Sample ON Uni(UO(()))) Unj Rx ~ Uni(SO(3)), rN ~ N(O, of,
Let A((rN, RN, N))
for n + N to do
Let l = n/N and Aot o (n/NN of((n - 1)/N) and similarly for Ao Ao R3, R33 y m s(xn, t);
Predict scores
Sample Zir, Zroos Zeor from (0, ((, (0,, Ao respectively;
Set Ar &- Aotar Zari
Set AR R(Do 2no)):
Set AO Ao + Zarr'
Compute Xn-1 f- A((Ar, AR, A0),
Return Xo:
C ARCHITECTURE DETAILS
As summarized in Section 4.5, we use convolutional networks based on tensor products of irre- ducible representations (irreps) of SO(3) [Thomas et al., 2018] as architecture for both the score and confidence models. In particular, these are implemented using the e3n library [Geiger et al., 2020]. Below, @w refers to the spherical tensor product of irreps with path weights w, and @ refers to normal vector addition (with possibly padded inputs). Features have multiple channels for each irrep. Both the architectures can be decomposed into three main parts: embedding layer, interaction layers, and output layer. We outline each of them below.
16

---

Published as a conference paper at ICLR 2023
C.I EMBEDDING LAYER
Geometric heterogeneous graph. Structures are represented as heterogeneous geometric graphs with nodes representing ligand (heavy) atoms, receptor residues (located in the position of the a- carbon atom), and receptor (heavy) atoms (only for the confidence model). Because of the high number of nodes involved, memory constraints. Moreover, sparsity can act as a useful inductive bias for the model, however, it is critical for the model to find the right pose that nodes that might have a strong interaction in the final pose to be connected during the diffusion process. Therefore, to build the radius graph, we connect nodes using cutoffs that are dependent on the types of nodes they are connecting:
necessary for the graph to sparsely connected for runtime and
Ligand atoms-ligand atoms, receptor atoms-receptor atoms, and ligand atoms-receptor atoms interactions all use cutoff of 5A, standard practice for atomic interactions. For the ligand atoms-ligand atoms interactions we also preserve the covalent bonds as sepa- rate edges with some initial embedding representing the bond type (single, double, triple and aromatic). For receptor atcms-receptor atoms interactions, we limit at 8 the maximum number of neighbors of each atom. Note that the ligand atoms-receptor atoms only appear in the confidence model where the final structure is already set.
2. Receptor residues-recepter residues use a cutoff of 15 Å with 24 as the maximum number neighbors for each residue.
3. Receptor paiduues-ligant atoms use a cutoff of 20 + 3 Oir À where Oer represents the current standard deviation of the diffusion translational noise present in each dimension (zero for the confidence model). Intuitively this guarantees that with high probability, any of the ligands and receptors that will be interacting in the final pose the diffusion model converges to are connected in the message passing at every step.
4. Finally, receptor residues are connected to the receptor atoms that form the corresponding amino-acid.
Node and edge featurization. For the receptor residues, we use the residue type feature as well as language model embedding obtained from ESM2 [Lin et al., 2022]. The ligand atoms have the following features: atomic number; chirality; degree; formal charge; implicit valence; the number of connected hydrogens; the number of radical electrons; hybridization type; whether or not it is in an aromatic ring; in how many rings is; and finally, 6 features for whether or not it is in a ring concatenated with sinusoidal embeddings of the diffusion time
of size 3, 4, 5, 6. 7, or 8. These (Vaswani et al., 2017] and, in the case of edges, radial basis embeddings of edge length [Schütt et al., 2017). These scalar features of each node and edge then transformed with learnable two-layer MLPs (different for each node and edge type) into a set of scalar features that are used as initial representations by the interaction layers.
Notation Let E), represent the heterogeneous graph, with V (Vi, V,) respectively ligand atoms and receptor residues (receptor atoms Va. present in the confidence model, are for simplicity not included here), and similarly & = (Eu,, Ere, Let h, be the node embeddings (initially only scalar channels) of node a, Cab the edge embeddings of (a, b), and (rras) radial basis embeddings of the edge length. Let of. a and of represent the variance of the diffusion kernel in each of the three components:
translational, rotational and torsional.
C.2 INTERACTION LAYERS
At each layer, for every pair of nodes in the graph, we construct messages using tensor products of the current node features with the spherical harmonic representations of the edge vector. The weights of this tensor product are computed based on the edge embeddings and the scalar features-denon the outgoing and incoming nodes. The messages are then aggregated at each node and used to update the current node features. For every node a of type ta:
h of of
h2 f- ha +
(t
Σ
((fat) @ea hg
(
rEfert
(16)
BEN?!
with whas ((ttart) (Eab, ho, hg)
17

---

Published as a conference paper at ICLR 2023
Here, indicates an arbitrary node type, Na {b (a,b) E Ei,,. the neighbors of a of type I, Y are the spherical harmonics up to = 2, and BN the (equivariant) batch normalisation. The orders of the output are restricted to a maximum of a dictionary of MLPs, which uses different sets of weights for different edge types (as an ordered pair so four types for the score model and nine for the confidence) and different rotational orders. Convolutional layers simultaneously operate with different sets of weights for different connection types (so sets of weights for the confidence model and 4 for the score at every layer) and generate scalar and vector representations for each node.
All learnable weights are contained in V,
C.3 OUTPUT LAYER
The ligand atom representations after the final interaction layer are used in the output layer to pro- duce the required outputs. This is where the score and confidence architecture differ significantly. On one hand, the score model's 's output is in the tangent space T,T @ Tx)0(3) @ (2)) This corresponds to having two SE(3) eqquivarian output vectors representing the translational and rotational score predictions and m S SE(3)-invaria) For each of these, design final tensor-product convolutions inspired by classical mechanics. On the other hand, the confidence model outputs single SE(3) (3n-iiiiaai scalar representing the con- fidence score. Below we detail how each of these outputs is generated.
output scalars representing the torsional score.
Translational and rotational scores. The translational and rotational score intuitively represent, respectively, the linear acceleration of the center of mass of the ligand and the angular acceleration of the rest of the molecule around the center. Considering the ligand as rigid object and given a set of forces and masses at each ligand, a tensor product convolution between the atoms and the center of mass would be capable of computing the desired quantities. Therefore, for each of the two outputs, we perform a convolution of each of the ligand atoms with the (unweighted) center of mass
1
Σ (frea VEVI
Sxea ha
[V
(17)
with thea ((((rca),)))) We restrict the output of v to a single odd and single even vectors (for each of the two scores). Since are using coarse- granied representations of the protein, the score will neither be even nor odd; therefore, we sum the even and odd vector representations of v. Finally, the magnitude (but not direction) of these vectors is adjusted with MIP taking as input the current magnitude and the sinusoidal embeddings of the diffusion time. Finally, we (revert the normalization) by multiplying the outputs by 1/o,, for the translational score and by the expected magnitude of a score in S0(3) with diffusion parameter roz (precomputed numerically).
Torsional score. To predict the m SE(3)-invarian scalar describing the torsional score, we use a pseudotorque layer similar to that of Jing et al. [2022]. This predicts scalar score 8r for each rotatable bond from the per-node outputs of the atomic convolution layers. For rotatable bond g (go.g) 91) and b E V/, let ligb and igg be the magnitude and direction of the vector connecting the center of bond g and b. We construct a convolutional filter T, for each bond g from the tensor product of the spherical harmonics with a 2 representation of the bond axis fg.
(18) @ is the full (i.e., unweighted) tensor product as described in Geiger & Smidt [2022], and the second term contains the spherical harmonics up to é 2 (as usual). This filter (which contains orders up to = 3) is then used to convolve with the representations of every neighbor on radius graph:
T (f) == (2 ( (f))
& = (ggbb) b) g rotatable bond, b E Vr
Cyb (()((((())))))
V(g.b) E E-
hg 1 Σ To((gh) Bra h W bEN,
(19)
with Ygb = T(egb, hg. h + hg)
Since the parity of the é = 2 spherical harmonic is even, this representation is indifferent to the choice of bond direction.
18

---

Published as a conference paper at ICLR 2023
bg E } and T( and I are MLPs with learnable parameters. Since unlike Jing et al. [2022], we use coarse-graines representations the parity also here neither even nor odd, the irreps in the output are restricted to arrays both even hy and odd hy scalars. Finally, produce a single scalar prediction for each bond:
Here, .Vg {b
5Tg = [1 + h
(20)
where a two-layer MLP with tanh nonlinearity and no biases. This is also "denormalized" by multiplying by the expected magnitude of score in SO(2) with diffusion parameter door
Confidence output. The single -in3n-invriiant scalar representing the confidence score output is instead obtained by concatenating the even and odd final scalar representation of each ligand atom, averaging these feature vectors among the different atoms, and finally applying a three layers MLP (with batch normalization).
D EXPERIMENTAL DETAILS
In general, all our code is available at https://gitti:: wttp://..oo/occoo////eee...///.. This in- cludes running the baselines, runtime calculations, training and inference scripts for DIFFDOCK, the PDB files of DIFFDock's predictions for all 363 complexes of the test set, and visualization videos of the reverse diffusion.
D.1 EXPERIMENTAL SETUP
Data. We use the molecular complexes in PDBBind [Liu et al., 2017] that were extracted from the Protein Data Bank (PDB) [Berman et al., 2003]. We employ the time-split of PDBBind proposed by Stärk et al. [2022] with 17k complexes from 2018 earlier for training/validation and 363 test struc- tures from 2019 with no ligand overlap with the training complexes. This is motivated by the further adoption of the same split [Lu et al., 2022] and the critical assessment of PDBBind splits by Volkov al. [2022] who favor temporal splits over artificial splits based on molecular scaffolds or protein sequence/struutuue similarity. For completeness, ilarity splits in Appendix F. We download the PDBBind data as it is provided by EquiBind from https: zenodo. org/ record/ee 6408497. These files before adding any potentially missing hydrogens, correcting hydrogens, and correctly flipping his- tidines with the reduce library available https : /gggttttt.///////// com r labuuke/ reduce.
also report the results protein sequence sim-
preprocessed with Open Babel
Metrics. To evaluate the generated complexes, we compute the heavy-atom RMSD between the predicted and the crystal ligand atoms when the protein structures are aligned. To account for per- mutation symmetries in the ligand, we use the symmetry-correcter RMSD of sPyRMSD [Meli & Biggin, 2020]. For these RMSD values, we report the percentage of predictions that have an RMSD that is less than 2A. We choose 2 since much prior work considers poses with an RMSD less that 2 as "good" or successful [Alhossary et al., 2015; Hassan al., 2017: McNutt et al., 2021]. This a chemically relevant metric, unlike the mean RMSD as detailed in Section since for further down- stream analyses such as determining function changes, a prediction is only useful below a certain RMSD error threshold. Less relevant metrics such as the mean RMSD are provided in Appendix F.
D.2 APO-STRUCTURE DOCKING
Although large and comprehensive, the PDBBind benchmark only evaluates the capacity that var- ious docking methods have bind ligands to their corresponding receptor colo-sttuccure. This is much simpler and less realistic scenario than what is typically encountered in real applications where docking for new ligands is done against apo or holo-structures bound to different ligand. In particular, since the development of accurate protein folding methods [Jumper et al., 2021]. docking programs are often run on top of Al-generated protein structures. With this in mind, we develop benchmark where combine the complex prediction of PDBBind with protein structures generated by ESMFold [Lin et al., 2022].
The structures obtained running esmfoldvv? the sequences extracted from the PDB pro- tein files of PDBBind. Protein complexes were passed to ESMFold by concatenating (with ' the sequences of the proteins and removing all water molecules and other ligands. Predictions were
19

---

Published as a conference paper at ICLR 2023
48GB A6000 GPUs, but for 12/361 complexes the prediction out of memory even when reducing the chunk.size and were, therefore, discarded.
The main design choice when generating this benchmark relies on how to best align the PDBBind complex with the ESMFold structure to obtain the "ground-truth docked prediction on the ESMFold structure. An unbiased global alignment of the two protein structures is not desirable because difference in structure not affecting the pocket where the ligand binds would cause the two pockets to misalign; on the other hand, only aligning residues within a single arbitrary pocket cutoff has many undesirable cases where too many or too few residues are selected or not weighted properly.
Instead, align receptors' residues with the Kabsch algorithm using exponential weighting. for every receptor x its weight is W'x = e x where À is smoothing factor and dx is the minimum distance of x to a ligand atom the original complex, this way residues closer to the ligand will have a higher weight in the alignment. For each complex, we individually select À E [0, 1] so that it preserves distances as best as possible, in partieular, we use the L-BFGS-B [Byrd et al., 1995] from scipy [Virtanen et al., 2020b] to minimize:
min
E XEX yEy
Σ
1
2
||x
||xe-Xl
[xx(A) yll
where ||x,- - y| and (Xx correspond to the distances between protein residue x and ligand atom y respectively in the original crystal structure from PDBBind and in the complex structure obtained aligning the ESMFold structure with smoothing parameter A. We use inverse distances to give more importance to residues closer to the ligand (in either structure) and avoid steric clashes. We only consider protein backbones because the side-chain predictions are often less reliable and their structure typically changes upon binding.
Thus we obtain protein structures on which we run the docking methods and the associated docked ligand positions that we use to evaluate them.
D.3 IMPLEMENTATION DETAILS: IIYPERPARAMETERS, TRAINING, AND RUNTIME MEASUREMENT
Training Details. We use Adam [Kingma & Ba, 2014] as optimizer for the diffusion and the con- fidence model. The diffusion model with which we run inference uses the exponential moving average of the weights during training, and we update the moving average after every optimization step with a decay factor of 0.999. The batch size is 16. We run inference with 20 denoising steps on 500 validation complexes every 5 epochs and use the set of weights with the highest percent- age of RMSDs less than 2 as the final diffusion model. We trained our final score model on four 48GB RTX A6000 GPUs for 850 epochs (around 18 days). The confidence model is trained on single 48GB GPU. For inference, only a single GPU is required. Scaling up the model size seems to improve performance and future work could explore whether this trend continues further. For the confidence model uses the validation cross cross-entropy loss is used for early stopping and training only takes 75 epochs. Code to reproduce all results including running the baselines to perform docking calculations for new complexes is available at https: / / //ggtttuut thub. com/ gcorso f fDock.
Hyperparameters. For determining the hyperparameters of DIFFDOCK's 's score model, we trained smaller models (3.97 million parameters) that fit into 48GB of GPU RAM before scaling it up to the final model (20.24 million parameters) that was trained on four 4SGB GPUs The smaller models were only trained for 250 or 300 epochs, and we used the fraction of predictions with an RMSD below 2 on the validation set choose the hyperparameters. Table 2 shows the main hyperparameters we tested and the final parameters of the large model we use to obtain our results. We only did little tuning for the minimum and maximum noise levels of the three components of the diffusion. For the translation, the maximum standard deviation is 19Å. We also experimented with second-order features for the Tensor Field Network but did not find them to help. The complete set of hyperparameters next to the main ones we describe here can be found in our repository. From the start we have divided the inference schedule into 20 time steps, the effect of using more or fewer steps for inference is discussed in Appendix F.3. As found that the large-scale diffusion models overfit the training data on low-levels of noise we stop the diffusion early after 18 steps. At the last diffusion step noise is added.
20

---

Published as a conference paper at ICLR 2023
The confidence model has 4.77 million parameters and the parameters tried are in Table 3. We generate 28 different training poses for the confidence model (for which it predicts whether or not they have RMSD below 2A) with small score model. The model used to generate the training samples for the confidence model does not need to be the same one that the model will be applied to at inference time.
Table 2: The hyperparameter options we searched through for DIFFDOCK's ' score model. This was done with small models before scaling up to a large model. The parameters shown here that impact model size (bottom half of the table) are those of the large model. The final parameters for the large DIFFDOCK model are marked bold.
PARAMETER
SEARCH SPACE
USING ATOMS THE PROTHIN GRAPH
YES, No
USING LANGUAGE MODEL. EMBEDDINGS
YES, No
USING LIGAND HYDROGENS
YES, No
USING EXPONENTIAL MOVING AVERAGE
YES, No
MAXIMUM NUMBER NEIGHBORS IN PROTEIN GRAPH
10. 16. 24. 30
MAXIMUM NEIGHBOR DISTANCE PROTEIN GRAPH
5. 10, 15, 18, 20.
DISTANCE EMBEDDING МETHOD
SINUSOIDAL, GAUSSIAN
DROPOUT
O. 0.05, 0.1, 0.2
LEARNING RATES
0.01, 0.008, 0.003, 0.001, 0.0008, 0.0001
BATCH SIZE
8. 16. 24 RELU
NON LINEARITIES
CONVOLUTION LAYERS
NUMBER SCALAR FEATURES
NUMBER OF VECTOR FEATURES
10
Table The hyperparameter options we searched through for DIFFDOCK's confidence model. The final parameters marked in bold.
PARAMETER
SEARCH SPACE
USING ALL ATOMS FOR THE PROTEIN GRAPH
YES, No YES, No No
EMBEDDINGS USING LIGAND HYDROGENS
LANGUAGE
USING EXPONENTIAL MOVING AVERAGE
No 10, 16, 24, 30
MAXIMUM NUMBER OF NEIGHBORS PROTEIN GRAPH
MAXIMUM NEIGHBOR DISTANCE PROTEIN GRAPH
5. 10. 15. 18, 20. 30 SINUSOIDAL 0, 0.05, 0.1, 0.2
DISTANCE EMBEDDING METHOD
DROPOUT
LEARNING RATES
0.03, 0.003, 0.0003, 0.00008 16
BATCH SIZE
LINEARITIES
RELU
CONVOLUTION LAYERS
5
NUMBER SCALAR FEATURES
24
NUMBER OF VECTOR FEATURES
Runtime. Similar to all the baselines, the preprocessing times are not included in the reported run- times. For DIFFDOCK the preprocessing time negligible compared the rest the inference time where multiple reverse diffusion steps are performed. Preprocessing mainly consists of forward pass of FSM2 to generate the protein language model embeddings, RDKit's conformer generation, and the conversion of the protein into radius graph. We measured the inference time when running on an RTX A100 40GB GPU when generating 10 samples. The runtimes we report for generating 40 samples and ranking them extrapolations where multiply the runtime for 10 samples by 4. In practice, this only gives an upper bound on the runtime with 40 samples, and the actual runtime should be faster.
Statistical significance. To determine the statistical significance the superior performance of our method we used the paired two-sample t-test implemented in sc py [Virtanen al., 2020a).
21

---

Published as a conference paper at ICLR 2023
D.4 BASELINES: IMPLEMENTATION USED SCRIPTS, AND RUNTIME DETAILS
Our scripts to run the baselines are available at ntttps://githa comm/goors//Diffooo For obtaining the runtimes of the different methods, we always used 16 CPUs except for GLIDE as explained below. The runtimes do not include any preprocessing time for any of the methods. For instance, the time that it takes run P2Rank is included for TANKBind, and P2Rank SMINA/GNINA since this receptor preparation only needs to be run once when docking many ligands to the same protein. In applications where different receptors are processed (such as reverse screening), the experienced runtimes for TANKBind and P2Rank SMINA/GNINA will thus be higher.
We note that for all these baselines we have used the default hyperparameters unless specified dif- ferently below. Modifyyng some of these hyperparameters (for example the scoring method's ex- haustiveness) will change the runtime and performance tradcoff's (e.g., if the searching routine is left running for longer then better poses are likely to be found), however, we leave these analyses to future work.
SMINA [Koes et al., 2013] improves Autodock Vina with a new scoring-function and user- friendliness. The default parameters were used with the exception of setting -unum modes 10. To define the search box, we use the automatic box creation option around the receptor with the default buffer of 4A all 6 sides.
GNINA [McNutt et al., 2021] builds on SMINA by additionally using a learned 3D CNN for scoring.
The default parameters were used with the exception of setting --nummmoeee mmmmmodes 10. To define the search box, we use the automatic box creation option around the receptor with the default buffer of 4A on all 6 sides.
Quick Vina-W [Hassan et al., 2017] extends the speed- optimized Quick Vina 2 [Alhossary et al., 2015] for blind docking. We reuse the numbers from Stärk al. [2022] which had used the default parameters except for increasing the exhaustiveness to 64. The files were preprocessed with the prepare 1 gand4 py and paulpres recommended by the Quick W authors.
py scripts of the MGLTools library is
Autodock Vina [Trott& & Olson, 2010] is older docking software that does not perform as well as the other more recent search-based baselines, but it is well-established tool. We reuse the numbers reported in TANKBind [Lu et al., 2022]
GLIDE [Halgren et al., 2004] is strong heavily used commercial docking tool. These methods all use biophysics based ecoring-functions We reuse the numbers from Stärk et al. [2022] since we do not have license. Running GLIDE involves running their command line tools for preprocessing the structures into the files required to run the docking algorithm. As explained by Stärk et al. [2022], the very high runtime of GLIDE with 1405 seconds per complex is partially explained by the fact that GLIDE only uses single thread when processing a complex. This fact and the parallelization options of GLIDE are explained here ht https: : / /www. sch rchrodi ngge... com/ ccom/bb//// 65. With GLIDE, it is possible to start data- datalparallel processes that compute the docking results for a different complex parallel. However, each process also requires a separate software license.
EquiBind [Stärk et al., 2022], we reuse the numbers reported in their paper and generate the predic- /ataarak///uiBa nd.
tions that we visualize with their code https : / //gi//ggtt github. com
TANKBind [Lu et al., 2022], we use the code associated with the paper at https: /ttp::////tttti tub. com/ Lwwellwwiii 10117/TaakBind The runtimes do not include the runtime of P2Rank or any pre- processing steps. In Table 1 we report two runtimes (0.72/2.5 sec). The first is the runtime when making only the top-1 prediction and the second for producing the top-5 predictions. Producing only the t predictions is faster since TANKBind produces distance predictions that need to be converted to coordinates with a gradient descent algorithm and this step only needs to be run once for the top-1 prediction, while it needs to be run 5 times for producing outputs. To obtain our runtimes we run the forward pass of TANKBind on GPU (0.28 seconds) with the default batch size of 5 that is used in their GitHub repository. To compute the time the distances-to co-ooosiio--ooinatteeee version step takes, we run the file baseline run tankbind parallel sh in our repository, Intel Xeon Gold
which parallelizes the computation across 16 processes which also
22

---

Published as a conference paper at ICLR 2023
6230 CPU. This way. obtain 0.44 seconds runtime for the conversion step of the top-1 prediction (averaged over the 363 complexes of the testset).
P2Rank [Krivák & Hoksza, 2018], is tool that predicts multiple binding pockets and ranks them. We use it for running TANKBind and P2Rank + SMINA/GNINA. We download the program from https: / //g thub. com/ rdk / p2rank and run it with its default parameters.
EquiBind + SMINA/GNINA [Stärk et al., 2022], the bounding box in which GNINA/SMINA searches for binding poses is constructed around the prediction of EquiBind with the --autobox__ 1 i gand option of GNINA/SMINA. EquiBind is thus used to find the binding pocket and SMINA/GNINA find the exact final binding pose. We --autobox__aatee additional 10À on all 6 sides of the bounding box following [Stärk et al., 2022].
10 to add an
EquiBind= SMINA /GNINA flez The same EquiBind + SMINA/GNINA but with the flex- ibility in terms of torsion angles activated for the sidechains that have an atom within 3.5A of the output from EquiBind.
P2Rank + SMINA/GNINA. The bounding box in which GNINA/SMINA searches for binding poses is constructed around the pocket center that P2Rank predicts as the most likely binding pocket. P2Rank is thus used to find the binding pocket and SMINA/GNINA to find the exact final binding pose. The diameter of the search box is the diameter of a ligand conformer generated by RDKit with an additional 10Å on all 6 sides of the bounding box.
P2Rank + SMINA fee IGNINA flez The same as P2Rank + SMINA/GNINA but with the flexibility in terms of torsion angles activated for the sidechains that have an atom within 3.5 the radius of the ligand away from the pocket center of P2Rank. Additionally we consider at most flexible sidechains with -- flex.max 10.
SMINA/GNINA + SMINA flour IGNINA flex The same as EquiBind + SMINA//NINA but the ini- tial ligand is provided by SMINA/GNINA ran without sidechain flexibility. The flexibility is in terms of torsion angles activated for the sidechains that have an atom within 3.5 of the output from SMINA/GNINA.
E ADDITIONAL DISCUSSION
E.1 ACCESS TO BOUND PROTEIN STRUCTURE
One limitation of DIFFDOCK is that it assumes access to the bound structure of the protein known as holo-structure. Although most of the literature in molecular docking makes this assumption, in practice, often only has access to the unbound apo protein structure the holo structure of the protein bound to a different ligand. Tackling the problem of binding to apo structures challenging due to the limited amount of data of unbound structures that have atom-to-atom correspondence to the holo structure with which they could be aligned. In DIFFDOCK the limitation of docking to holo structures may be less pronounced than for the search based docking methods that we compare against since DIFFDOCK'S 's score model only uses the positions of the alpha carbon atoms (and not the side chain atoms). Thus DIFFDOCK would also work well for binding to apo structures when most of the conformational change during binding lies in the side chains and the backbone stays mostly rigid. However, in order to fully model binding to apo structures, one needs to additionally model protein flexibility, we leave the full treatment of this problem to future work.
E.2 TORSIONAL VS EUCLIDEAN FLEXIBILITY
A ligand pose (and more generally speaking, a molecule conformation) consists of position in 3D-space 3 for each of n atoms of the molecule and can thus be considered as an element of R33 Each such pose (or conformation), can alternatively be described in terms of its bond lengths, angles, and torsion angles (see Jing et al. [2022], Appendix for a more formal discussion). Be- cause of the nature of covalent interactions, bond lengths and angles (local structures) are highly energetically constrained, and any specific bond length or angle in a molecule will take on only a very narrow range of values, whether in isolation or bound to a protein receptor. Thus, ligand poses or conformations in which bond lengths or angles are strained differ significantly from their standard values-cal can be easily judged to be energetically unfavorable. On the other hand, the en-
23

---

Published as a conference paper at ICLR 2023
ergetic profile of varying torsion angles is significantly smoother, they depend in large part weaker, noncovalent interactions. A change chemical environment, or an interaction with another molecule (such protein) can alter these profiles. Thus, say that there is significant flexibility in the torsion angles, in contrast to the more rigid local structures.
Since bond lengths and angles are highly constrained, the set of chemically plausible ligand poses a very small subset of all possible assignments of atoms to R that is, a very restricted subset of 3 The set of ligand poses that satisfy all the bond length and angle constraints with equality appears as a lower-dimensional manifold. To a good approximation. all actual ligand poses lie on this manifold. Thus, we are interested in developing generative model on this manifold, as any sample outside of it is chemically nonsensical.
Futher, the constrained values of standard bond lengths and angles are easily and quickly predicted using standard cheminformatics routines, such as those in RDKit. Indeed, far simpler rules of thumb fixed length for each type of bonding pair, and fixed angles determined solely by molecu- lar geometry-art already quite accurate and widely taught. Consequently, the manifold of plausible poses is very easy to find-one merely has to generate any plausible conformation or ligand pose with RDKit. In previous work on conformer generation, Jing et al. [2022] verified that the manifolds corresponding to local structures generated by RDKit are, on average, less than 0.5 À RMSD away from the true conformations.
The facts that (1) the space of plausible ligand poses is described by a manifold and (2) this manifold is easy to find, motivate the development of diffusion generative model the manifold rather than the full ambient space R This significantly reduces the dimensionality of the state space, and thus of the function that the score network needs to approximate.
F ADDITIONAL RESULTS
F.1 PHYSICALLY PLAUSIBLE PREDICTIONS
Table 4: Steric clashes. Percentage of test complexes for which the predictions of the different methods exhibit steric clashes. Search-based methods never produced steric clashes.
Top-1
Top-5
Method
% steric clashes
% steric clashes
EQUIBIND
26
TANK BIND
6.6
3.6
DIFF (10)
2.8
DIFFDOCK (40)
2.2
2.2
Due to the averaging phenomenon of regression-based methods such TANKBind and EquiBind, they make predictions at the mean of the distribution. If alcatoric uncertainty is present, such as in case of symmetric complexes, this leads predicting the ligand be an un-physical state in the middle of the possible binding pockets as visualized in Figure 13, The Figure also illustrates how DIFFDOCK does not suffer from this issue and is able to accurately sample from the modes.
In the scenario when epistemic uncertainty about the correct ligand conformation is present, this often results in "squashed-up" predictions of the regression based methods as visualized in Figure 4. If there is uncertainty about the correct conformer, the square enror minimizing option is to put all atoms close the mean.
These averaging phenomena in the presence of either aleatoric or epistemic uncertainty cause the regression-based methods to often generate steric clashes and self intersections. To investigate this quantitatively, we determine the fraction of test complexes for which the methods exhibit steric clashes. We define a ligand as exhibiting a steric clash if one of its heavy atoms within 0.4À of heavy receptor atom. This cutoff is used by protein quality assessment tools and in previous literature [Ramachandran et al., 2011]. Table shows that DIFFDOCK, as generative model, pro- duces fewer steric clashes than the regression-based baselines. We generally observe no unphysical
24

---

Published as a conference paper at ICLR 2023
Ligand self-intersections. TANKBind (blue), EquiBind (cyan), DIFFDOCK (red), and crystal structure (green). Due the averaging phenomenon that occurs when epistemic uncertainty present, the regression-based deep learning models tend to produce ligands with atoms that are close together, leading to self-intersections. DIFFDOCK, as generative model, does not suffer from this averaging phenomenon, and never found self-intersectii in any of the investigated results of DIFEDOCK.
Figure 4:
Figure 5: Chemically plausible local structures. TANKBind (bluc), EquiBind (cyan), and DIFF- DOCK (red) structures for complex 6g EquiBind (without their correction step) produces very unrealistic local structures and TANKBind, e.g., produces non non-plane planar aromatic rings. DIFFDOcK's local structures are the realistic local structures of RDKit.
predictions from DIFFDOCK unlike the self intersections that, e.g., TANKBind produces (Figure 4) or its incorrect local structures (Figure 5). This is also visible in the randomly chosen examples of Figure 12 and can be examined in our repository, where we provide all predictions of DIFFDOCK for the test set.
25

---

Published as a conference paper at ICLR 2023
F.2 FURTHER RESULTS AND METRICS
In this section, we present further evaluation metrics on the results presented in Table 1. In particular, for both top-1 (Table 5) and top-5 (Table 6) we report: 25th, 50th and 75th percentiles, the propor- tion below 2 and below 5 of both ligand RMSD and centroid distance. Moreover, while Volkov al. [2022] advocated against artificial protein set splits and for time-based splits, for complete- ness, in Table 7 and Figure 7, we report the performances of the different methods when evaluated exclusively on the portion of the test set where the UniProt IDs of the proteins are not contained in the data that is seen by DIFFDOCK in its training and validation.
To assess the impact of the molecule size on the performance of DIFFDOCK and GNINA, we provide scatter plots in Figure 8 showing that the correlation between RMSD and the number of rotatable bonds or the number of atoms the molecule is similar for both methods. Figure we show that, based Tanimoto similarity, the performance of DIFFDOCK does not depend the test ligand's similarity to the ligands that have already been seen during training. The Spearman rank correlation coefficient between the RMSD and the Tanimoto similarity to the closest ligand in the training set is negligible with being -0.031.
Further, in Table setting when the flexibility of sidechains in the pocket is turned on. These show that simply adding sidechain flexibility does help these methods. Finally, in Figure 10, we plot the relationship between the quality of the ESMFold structure and the performance and holo-structures. While is approximately correct, very small variations in the apo-structure backbone (and sidechains) causes search-based methods like GNINA almost never find the right pose.
provide the performance of SMINA and GNINA in the apo-structure docking
the docking methods on apo retains large part of its accuracy when the protein backbone
Fraction with lower RMSD
0.7
DifiDcck
Fraction with lower RMSD
0.7
Method
0.6
GNINA
GUDDE
0.6
DiffDock
GNINA
0.5
SMINA
0.5
SMINA
0.4
TANKBind
EquiBind TANKBind
EquiBind
0.3
0.3
0.2
0.2
RMSD (A)
RMSD (A)
Figure 6: Cumulative density histogram of the methods' RMSD: left on holo crystal structures, right on apo ESMFold structures.
Fraction with lower RMSD
DiffDock
Fraction with RMSD < -2
1
Baseline performance
GLIDE
Perfect selection
0.6
GNINA
Confidence
SMINA
OVinaW
0.6
0.4
Tand
EquiBind
0.3
0.4
0.2
0.1
40
60
RMSD (A)
Percentage of rejected complexes
Figure 7: PDBBind docking unseen receptors. Left: cumulative density histogram of the meth- ods' RMSD. Right: Percentage of predictions with RMSD below 2Á when only making predictions for the portion of the dataset where DIFFDOCK is most confident.
26

---

Published as a conference paper at ICLR 2023
Table 5: Top-1 PDBBind docking.
Ligand RMSD
Centroid Distance
% below threshold f
% below thresh.
Percentiles
Percentiles
Methods
25th
50th
75th
5 5
A
25th
50th
75th
5 Á
2
AUTODOCK VINA
5.7
10.7
21.4
21.2
5.5
1.9
6.2
20.1
47.1
26.5
QVINA-W
2.5
7.7
23.7
40.2
20.9
0.9
3.7
22.9
54.6
41.0
GNINA
2.4
17.9
40.8
53.6
SMINA
17.9
38.0
18.7
1.0
59.8
41.6
GLIDE (c.)
2.6
9.3
28.1
33.6
21.8
0.8
5.6
26,9
48.7
36.1
EQUIBIND
3.8
6.2
10.3
39.1
5.5
1.3
2.6
7.4
67.5
40.0
TANKBIND
2.5
4.0
8.5
59.0
20.4
0.9
1.8
4.4
77.1
55.1
PIRANK+SMINA
2.9
6.9
16.0
43.0
20.4
0.8
2.6
14.8
60.1
44.1
GNNAA
1.7
5.5
15.9
47.8
28.8
0.6
2.2
14.6
60.9
48.3
EQUIBIND+SMINA
11.2
43.6
EQUIBIND+GNINA
1.8
4.9
13
50.3
28.8
0.6
1.9
9.9
66.5
50.8
DICK (10)
61.7
35.0
80.7
63.1
DOCK (40)
1.4
3.3
7.3
63.2
38.2
0.5
1.2
3.2
80.5
64.5
Table 6: Top-5 PDBBind docking.
Ligand RMSD
Centroid Distance
% below threshold f
% below thresh.
Percentiles
Percentiles
Methods
25th
50th
75th
5 Á
2 A
25th
50th
75th
5 5
2
GNINA
1.6
11.8
52.8
29.3
0.6
2.0
8.2
66.8
49.7
SMINA
1.7
4.6
9.7
53.1
29.3
0.6
1.85
6.2
72.9
50.8
TANKBIND
2.1
3.4
67.5
24.5
0.8
1.4
2.9
86.8
62.0
P2RANK+SMINA
1.5
4,4
14.1
54.8
33.2
0.6
1.8
12.3
66.2
53.4
GNNINA
1.4
3.4
12.5
60.3
38.3
0.5
1.4
9.2
69.3
57.3
EQUIBIND+SMINA
1.3
3.4
8.1
60.6
38.6
0.5
1.3
5.1
74.9
58.9
EQUIBIND+GNINA
1.4
61.7
73.7
DIFFDOCK (10)
1.2
2.7
4.9
75.1
40.7
0.5
1.0
2.2
87.0
72.3
DIFKDOCK (40)
1.2
2.4
5.0
44.7
0.4
0.9
1.9
88.0
76.7
Table PDBBind docking on unseen receptors. Percentage of predictions for which the RMSD to the crystal structure is below 2 and the median RMSD. indicates the method run exclusively on CPU, "_" means not applicable; some cells are empty due infrastructure constraints.
Top- %<2
RMSD
Top-5
RMSD
Average Runtime (s)
Method
Med.
%<2
Med.
AUTODOCK VINA
1.4
16.6
205*
QVINAW
15.3
10.3
49#
GNINA
14.0
13.6
23.0
7.0
127
SMINA
14.0
8.5
21.7
6.7
126*
GLIDE
19.6
18.0
1405*
EQUIBIND
0.7
9.1
0.04
TANKBIND
6.3
5.0
11.1
4.4
0.7/2.5
DIFFDOCK (10)
15.7
6.1
21.8
4.2
10
DIFFDOCK (40)
20.8
6.2
28.7
3.9
40
27

---

Published as a conference paper at ICLR 2023
RMSD crystal structure
RMSD to crystal structure
4-
DiffDock
DifiDock
20
40
60
80
10
20
30
40
50
60
Number heavy atoms
Number of rotatable bonds (without
Figure 8: Relationship between the number of atoms and RMSD. Left: Seatter plot of the RMSD and the number of atoms. Right: Scatter plot of the RMSD and the number of rotatable bonds in the ligand.
10
RMSD to crystal structure
8-
4
2
50
100
150
200
Rank by Maximum Tanimoto Similarity
Figure 9: Scatterplot for the RMSD and the Tanimoto similarity of the test ligand to the closest ligand in the training data, The order by maximum Tanimoto similarity such that the lowest similarity is at the left and the highest similarity is at the right on the x-axis.
Table 8: PDBBind blind docking on apo proteins with sidechain flexibility baselines. All methods receive small molecule and are tasked to find its binding location, orientation, and formation. Shown is the percentage predictions with RMSD < 2Á and the median RMSD (see Appendix D.3). In parenthesis, we specify the number of poses sampled from the generative model. indicates that the method runs exclusively on CPU. No method received further training or tuning on ESMFold structures.
Apo ESMFold proteins
Top-1
RMSD
Top-5 RMSD
Average Runtime (s)
Method
%x2
Med.
%x2
Med.
P2RANK+SMINA flaz
5.7
9.6
13.4
6.3
292*
P2RANK+GNINA flezz
8.3
10.9
13.4
6.4
294
EQUIBIND+SMINA flez
4.3
7.3
11.7
5.8
1145*
EQUIBIND+GNINA flez
6.6
9.8
14.6
6.1
1208
SMINA+SMINA fice
3.4
12.6
8.3
11.6
1145*
GNINA+GNINA fire
1.7
22.1
5.1
20.0
1208
DIFF (10)
21.7
5.0
31.9
3.3
DOCK (40)
20.3
5.1
31.3
3.3
40
28

---

Published as a conference paper at ICLR 2023
DifiDock PDB DiffDock ESMFold
Fraction with RMSD < 2A
GNINA ESMFold
H
< 0.5 A
A < 1.5 Backbone pocket RMSD
1.5 A
Figure 10: Relationship between the ESMFold and docking accuracy. The test set complexes on PDBBind are divided into three similar-sized groups based on the RMSD of the pocket residues (defined as i 10 À to ligand atom) in the aligned ESMFold generated structure. The performance of DIFFDOCK and GNINA are shown when docking to both the ground truth crystal structure and the ESMFold structure. The baseline methods' performance on crystal structures is also negatively correlated with the ESMFold accuracy because the complexes where the methods do badly tend to be larger and with fewer examples in PDB(Bind).
F.3 ABLATION STUDIES
Below we report the performance of our method over different hyperparameter settings. In partic- ular, we highlight the different ways in which it is possible to control the tradeoff between runtime and accuracy in our method. These mainly are: (I) model size, (2) diffusion time, and (3) diffusion samples.
Model size. The final DIFFDOCK score model has 20.24 million parameters from its convolution layers with 48 scalar and 10 vector features. In Table 9 show the results for smaller score model with 5 convolutions, 24 scalar, and 6 vector features resulting in 3.97 million parameters that can be trained on single 48GB GPU. The confidence model used is the same for both score models. We find that scaling up the model size helped improve performance which we did as far as possible using four 48GB GPUs for training. Scaling the model size further is a promising avenue for future work.
Protein embeddings. As described in Section 4.5 and Appendix C, the architecture features of protein residues the language model embeddings from ESM2 [Lin et al., 2022] in order for the model more easily reason about the protein sequence. In Table we show that while these provide some improvements they are not necessary to obtain statc-of-the-arr performance.
initial
Diffusion steps. Another hyperparameter determining the runtime of the method during inference is the number of steps we take during the reverse diffusion. Since these are applied sequentially DIFFDOCK's runtime scales approximately linearly with the number of diffusion steps. the rest of the paper, we always use 20 steps, but in Figure 11 we show how the performance of the model varies with the number of steps. We note that the model reaches nearly the full performance even with just 10 steps, suggesting that the model can be sped up 2x with small drop in accuracy.
Diffusion samples. Given a score-based model and a number of steps for the diffusion model, remains to be determined how many independent samples N to query from the diffusion model and then feed to the confidence model. As expected the more samples the confidence model receives the more likely it is that it will find pose that it is confident about and, therefore, the higher the performance. The runtime of DIFF DOCK on GPU scales sublinearly until the different samples fit in parallel in the model (depends the protein size and the GPU memory) and approximately linearly
29

---

Published as a conference paper at ICLR 2023
Table Model size and protein embeddings comparison. All methods receive small molecule and are tasked to find its binding location, orientation, and conformation. Shown is the percentage of predictions for which the RMSD to the crystal structure below 2 and the median RMSD.
Top-1 %<2
RMSD (A)
Top-5 %<
RMSD (A)
Average
Method
Med.
Med.
Runtime (s)
DIFFDOCK-SMMALLNOOESS (10)
26.2
4.7
32.0
3.2
DIFDOCKK-MMALLNNOSSS (40)
28.4
3.8
37.7
2.6
28
DIFEDOCK-SMALL (10)
26.0
4.3
33,3
3.2
7
DIFFDOCK-SMALL (40)
31.1
4.0
38.0
2.7
28
DIFFDOCK-NOESM (10)
33.9
3.8
39.4
2.8
10
DIFDDOCK-NOESM (40)
34.2
3.5
42.7
2.4
40
DIFFDOCK (10)
35.0
3.6
40.7
2.7
10
DIFFDOCK (40)
38.2
3.3
44.7
2.4
40
0.5
Fraction with RMSD 2A
0.4
0.3
0.2
Top-1 performance
0.1
Top-5 performance
Top-10 performance
10
20
30
40
50
60
Number of reverse diffusion steps
Figure 11: Ablation study on the number of reverse diffusion steps.
for larger sample sizes (however it can be easily parallelized across different GPUs). In Figure we show how the success rate for the top-1, top-5, and top- 10 prediction change as a function N. For example, for the top-1 prediction, the proportion of the prediction with RMSD below 2A varies between 22% of a random sample of the diffusion model (N = 1) to 38% when the confidence model is allowed to choose between 40 samples.
30

---

Published as a conference paper at ICLR 2023
F.4 VISUALIZATIONS
Randomly picked examples. The predictions of TANKBind (blue), EquiBind (cyan), GNINA (magenta), DIFFDOCK (red), and crystal structure (green). Shown are the predictions once with the protein and without below. The complexes were chosen with random number generator from the test set. TANKBind often produces self intersections (examples at the top-right; middle- middle; middle-right; bottom-right). DIFFDOCK and GNINA sometimes almost perfectly predict the bound structure (e.g., top-middle). The complexes in reading order are: 6p8y, 6mo8, 6pya, 6t6a, 6e30, 6hld, 6qzh, 6hhg. 6qln.
Figure 12:
31

---

Published as a conference paper at ICLR 2023
Figure 13: Symmetric complexes and multiple modes. EquiBind (cyan), DIFFDOCK highest confidence sample (red), all other DIFFDOCK samples (orange), and the crystal structure (green). We see that, since is generative model. DIFFDOCK able produce multiple correct modes and to sample around them. Meanwhile, as a regression-baut based model, EquiBind is only able to predict a structure at the mean the modes. The complexes are unseen during training. The PDB IDs in reading order: 6agt, 6gdy, 6ckl, 6dz3.
32

---

Published as a conference paper at ICLR 2023
Figure 14: Reverse Diffusion. Reverse diffusion of a randomly picked complex from the test set. Shown are DIFFDOCK highest confidence sample (red), all other DIFFDOCK samples (orange), and the crystal structure (green). Shown are the 20 steps of the reverse diffusion process (in reading order) of DIFFDOCK for the complex 6oxx. Videos of the reverse diffusion are available at https: //github. http:////cooooooooooooo///////aa zanions/ README, mdd
33


Reviews:
[{"review_id": "YhbLnS_1_sn", "reviewer": "Reviewer_5meX", "summary": "", "questions": "", "limitations": "", "rating": 10, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "DiffDock is a diffusion model for molecular docking. The authors of the paper formulate the docking problem as a generative modeling task. In order to preserve the key conformational constraints of the molecules, compounds are mapped to a product space of the degrees of freedom, i.e. translation, rotation, and torsional angles. In such created space a diffusion procedure can be defined after ensuring that torsions are orthogonal to the roto-translations, and the mapping is bijective. The diffusion is defined in the product space with a stochastic differential equation and a score model. Furthermore, a confidence model is trained using RMSD values of docked compounds. The experiments show advantage of the proposed method over other classical and neural dockers in the blind docking task. ", "summary_of_the_review": "Based on the above comments, I recommend the acceptance of this paper.", "strength_and_weaknesses": "Strengths:\n- The authors indicate the problems with regression-based models and explain the motivation of using diffusion models in the context of pose prediction.\n- The problem definition and the execution of the model are very good. The authors notice that molecular conformations lie in a submanifold where all the operations should be performed instead of moving the atom positions freely.\n- The torsion updates are carefully defined so that they possess properties described in Propositions 1 and 2.\n- A confidence model is trained as a classification model that predicts whether the RMSD of the generated pose is less than 2A.\n- SE(3)-invariant models are used for the score model and the confidence model. This is important to distinguish two reflections of the same compound (enantiomers).\n- The experimental results show the advantages of the proposed method in terms of binding pose prediction and inference runtime.\n- Additional experiments as well as all the training details are provided in the supplementary materials.\n\nWeaknesses:\n- The only problem I see in the definition of the product space is that it still allows for the steric clashes (either with the protein or other ligand atoms during the torsional rotations). In Figure 5, you argue that self-intersections are not observed for your model. Do you have any intuition why these phenomena do not occur for DiffDock?\n- If the confidence model is trained in a classification setup, then the output probabilities do not account for the uncertainty of this model. Do you think it could break the ranking for the unseen binding pockets and ligand poses? Could you elaborate on why you decided to use classification with the threshold of 2A?\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "Tw1-EuIvtBZ", "reviewer": "Reviewer_Yi6P", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors present a molecular docking method based on diffusion model. They suggest a diffusion process for this space and an architecture that learns on this process. The authors present results on a version of PDBbind. ", "summary_of_the_review": "The novelty of the proposed method is low because similar works using diffusion processes for molecules have been suggested before. The relevance for machine learning is very limited because many settings of the application domain have been used, therefore, it is unclear how to transfer this to similar ML problems. Technically, there are many severe problems, such as the limited experiments, problems in the method comparison, inconsistencies with previous method comparisons, missing error bars, and sloppy mathematical formalisms. Together with the many false and tenuous claims, the work defies some basic scientific practices.\n\nHowever, with some improvements this work will be highly interesting for the molecular docking community and will be a valuable contribution there.\n", "strength_and_weaknesses": "Strengths:\n\n* Clearly written and structured; hyperparameter selection clearly described\n* The principle idea to use diffusion models to predict 3D coordinates of atoms of a molecule is good (however not new; see weaknesses)\n\nWeaknesses\n\n* the work is not novel, but the use of diffusions for molecule conformations has been suggested before (see ref [1], [7])\n* the work is poorly embedded into related work and almost fully ignores decades of work in molecular docking\n* the mathematical formulation is sloppy and potentially wrong (see below), with notation overloads and inconsistencies\n* there is only a single experiment on one dataset which limits the conclusions\n* the work is strongly focused on a specific problem in an application domain, with no idea put forward how this would be relevant for ML in general\n* the method comparison neglects scientific standards to perform repetitions, provide error bars or confidence intervals and statistical tests. Almost all performance metrics are shown without error bars\n* the paper contains many ad-hoc and un-justified decisions\n* the paper contains many incorrect or tenuous claims\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "CzREHZd4RUv", "reviewer": "Reviewer_UKAX", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors propose DiffDock, an approach to global molecular docking as a generative modeling problem that aims to map the non-Euclidean manifold of ligand poses to the degrees of freedom involved in docking sampling by developing a diffusion process on this space. They propose an objective formalism that aims to maximize the expected proportion of predictions with RMSD < some tolerance ϵ motivating the training of a generative model that learns a distribution over ligand poses conditioned on the protein structure as context. In addition, a confidence model is proposed which is trained on the poses sampled by the diffusion model to rank them based on its confidence that they are within an error tolerance. Therefore posing and scoring is achieved by this process. The paper evaluates PDBBind dataset with the task of blind docking evaluated by the RMSD compared to the reference co-crystal structure. DiffDock is then compared to SOTA search-based and deep-learning methods. When evaluating the top scoring pose being <2A˚ from the co-crystal ligand, DiffDock **significantly** outperforms all compared methods with also being 3 to 12 times faster than the best search-based method.", "summary_of_the_review": "The authors propose DiffDock, a generative diffusion model tailored to the task of molecular docking. The approach fundamentally differs from regression-based frameworks and leverages a generative modeling approach to capture a diffusion process over the manifold describing ligand pose transformations and docking degrees of freedom. DiffDock demonstrates significant improvements over SOTA search-based and deep-learning approaches to global docking on the PDBBind dataset. It is not clear how exactly the protein-ligand structure data was curated and prepared from the PDBBind dataset or why the PDBBind database was chosen for the task of global docking. Still, the data was presented clearly, and the theory was well described. As outlined above, there are clear gaps that must be addressed before DiffDock will be useful for real-world applicability in drug discovery, and it is clear that there is a lack of domain expertise in computational chemistry and structure-based drug design. However, I believe the theoretical contributions and novelty of the manuscript warrant acceptance to ICML following revision.", "strength_and_weaknesses": "**Strengths**:\n\nThe proposed diffusion generative model formalism and application to the problem solving a global docking objective is novel and theoretically sound. The addition of the score-based generative model as a ranking mechanism takes the applicability of the diffusion probabilistic method to the next level in terms of posing, scoring, and creating a useful tool.\n\nThe comparison to search-based methods (as a quasi-null model) exemplified the significance of the results presented in the manuscript, and clearly DiffDock provides a significant improvement with the experimental protocol performed by the authors. The plots showing the fraction of RMSD vs. other factors further demonstrate the significant improvements compared to other methods and how well the distributions capture the data-generating process. \n\nWeaknesses:\nA large weakness of the work appears to be the lack of curation and preparation of the ligand and protein structures obtained from PDBBind. Structure-based drug design is only useful if the structures have been properly prepared by quantifying ionization/tautomeric states for ligands and H-bond optimization/protonation for proteins (garbage in~garbage out). When evaluating the code, it was clear that raw structures were input with minimal curation/preparation. If these raw structures were used as input to the other docking methods, then it is no surprise they performed poorly. The authors should describe their protein-ligand structure preparation protocol in detail. \n\nIt is not clear if PDBBind is the best dataset to use for evaluating the global docking problem, as all structures are holo-structures (co-crystal structures; i.e. structures in the ligand-bound conformation). A more relevant test set would be an apo (no ligand bound) structure of the same protein, similar to what is available in the PSCDB (http://idp1.force.cs.is.nagoya-u.ac.jp/pscdb). While this task would be more challenging, it would be more relevant for real-world problems encountered in drug discovery. It is well known in the docking community that predicting <2A˚ poses to holo-structures is much easier if the protein and ligand are prepared correctly, but recapitulating the crystallographic pose when docking to apo structures is much more challenging, even when integrating orthogonal biophysical data collected experimentally.\n\nIt appears that only Cα atoms are used for the node and edge featurization of the protein. This is a large limitation as structure-based drug design requires explicit representation of rotamers for amino acids as the orientation/protonation of these residues can dramatically alter what interactions can be made and ultimately whether a molecule will bind with sufficient potency to warrant detection in an in vitro assay.\n\nSome of the molecules generated by this approach are not physically possible to produce and would be incredibly unstable on planet earth.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "00j1bYZYT1", "reviewer": "Reviewer_pHbo", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper proposes a new way for molecule docking, which is to use diffusion model to form the docking pose prediction as a generation problem. In this work, the prediction is mainly based on the 3D translation group, rotation group and torsion angle change. The authors first critique the regression-based ligand propose prediction. Then define the ligand pose transformations, and the diffusion process, and then the training, model arch. The experiments on PDBBind dataset show that the method is strong with good verification results. The docking accuracy is highly improved and the running time is improved compared with some other methods. \n\nOverall speaking, this paper is strong in terms of performance, but I have lots of concerns.", "summary_of_the_review": "N/A", "strength_and_weaknesses": "Strengths:\n1. The authors think of the pose prediction task as a generation task and use the advanced diffusion model to perform the task definition. They also talked about the potential problem of the regression-based method. The regression-based method can not fit the overall multiple correct poses in a pocket, while the generative-based model can handle and provide different samples. A confidence model is used to keep the choice of the best ones of predictions.\n2. The authors consider the rotation, translation, and torsion angle change as the key messages that are variant during the generation, which largely reduces the freedom space of the ligand pose generation. \n3. The experiments are strong and effective compared to previous works, both in accuracy and running time. \n\nWeaknesses:\n1. Overall speaking, this paper is not very easy to read and follow (at least for me). The organization is good but not the details. The authors are experts in this domain (this paper is already put on arxiv and the authors are known), but the writing in this paper contains lots of parts that may not be friendly to the general readers. I can understand that there are lots of constraints or domain knowledge that the authors know, but for a paper, it should be easy for the readers to understand. Why I feel hard is that lots of items (words) and knowledge should be discovered by the ML readers. For a reviewer like me (I should say that I know basic knowledge about docking but am not so familiar, at least not familiar as authors), it is really hard for me to capture all of the detailed things that the authors are saying. In comparison, I read TankBand, which is much easier to follow and understand. I have also read the previous work Equiband. I strongly encourage the authors can use much simpler descriptions of the method. By the way, a simple and effective method is acknowledged in the ML committee I think. \n2. For some details, there lacks of explanations at all. For example, the m+6 degrees, 6 refers to the rototranslations, what is the detail of this 6? Product space, group space, aleatoric uncertainty, epistemic uncertainty, authors look like to define lots of new items. I am not sure whether these are common knowledge or defined by the authors only. I do feel that these make the readers hard to follow. Even for the diffusion method, details are missing. It's hard to say that 4.3 contains the necessary information for diffusion generation. Lots of details are put in the appendix. This is limited by the page constraints of ICLR. But if this is the reason, a journal that is without page length constraints is better? In the current main text, the training and inference is also not clear, which are placed in the appendix. As for section 3, I am not sure whether this really impresses me, the authors put the whole section to criticize the regression method. What I understand best is that the regression method can not handle the multiple correct poses problem, and the prediction is only the mean pose of these poses. The variance is not enough compared to the generation method. In my view, all of these can be said about the difference between the general regression model and the diffusion model. VAE, flow, diffusion, these models are all good at variance prediction. Is that necessary to put the whole section here, and the cases in figure 2 are somehow cherry-picked? \n3. In appendix, the training and inference procedure, the authors talk much about the problem of the different gaps caused by the initial conformation c. In my view, the description can be shortened since contains much about intuitively speaking, and conceptually speaking. Too much description also confused me instead of letting me understand better (only in my view when I read this part).\n4. When I read the details of the implementation, I have more concerns. The protein model is initialized of ESM2, a strong protein pre-trained model, but the authors are not mentioned in the main text at all. This strongly pre-trained model may help a lot. Besides, the module ending is also different, the authors include a lot of knowledge features of the molecule, for example, the number of rings and the details of the rings. there are all different from previous works. I am not sure whether these are the main contribution to the accuracy or the diffusion model. I would like to see fair comparisons of the model and encoding ways. \n5. For the small model and the large model, the small model is only used for the hyperparameter search? Even like this, it seems that the training time is also long, 18 days on 4 A6000 cards. And, why the inference is conducted on an A100 card then? ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Accept: poster