Equivariant Scalar Fields for Molecular Docking with Fast Fourier Transforms
Authors: Bowen Jing
Tommi S. Jaakkola
Bonnie Berger
Venue: ICLR 2024 poster (2024)
Keywords: protein structure
structural biology
drug discovery
molecular docking
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)

Abstract: Molecular docking is critical to structure-based virtual screening, yet the throughput of such workflows is limited by the expensive optimization of scoring functions involved in most docking algorithms. We explore how machine learning can accelerate this process by learning a scoring function with a functional form that allows for more rapid optimization. Specifically, we define the scoring function to be the cross-correlation of multi-channel ligand and protein scalar fields parameterized by equivariant graph neural networks, enabling rapid optimization over rigid-body degrees of freedom with fast Fourier transforms. The runtime of our approach can be amortized at several levels of abstraction, and is particularly favorable for virtual screening settings with a common binding pocket. We benchmark our scoring functions on two simplified docking-related tasks: decoy pose scoring and rigid conformer docking. Our method attains similar but faster performance on crystal structures compared to the widely-used Vina and Gnina scoring functions, and is more robust on computationally predicted structures. Code is available at https://github.com/bjing2016/scalar-fields.

Published as a conference paper at ICLR 2024
EQUIVARIANT SCALAR FIELDS FOR MOLECULAR DOCKING WITH FAST FOURIER TRANSFORMS
Bowen Jing', Tommi Jaakkola', Bonnie Berger's
CCSALL, Massachusetts Institute of Technology
Dept. of Mathematics, Massachusetts Institute of Technology bj Lingemit.edun {tommi, bab pelecsail.mit...mm
ABSTRACT
Molecular docking is critical to structure-based virtual screening, yet the through- put of such workflows limited by the expensive optimization of scoring func- tions involved in most docking algorithms. We explore how machine learning can accelerate this process by learning a scoring function with a functional form that allows for more rapid optimization. Specifically, define the scoring func- tion to be the cross-cores coss-correlltion of multi-channel ligand and protein scalar fields parameterized by equivariant graph neural networks, enabling rapid optimization over rigid-body degrees of freedom with fast Fourier transforms. The runtime of our approach can be amortized at several levels of abstraction, and is particularly favorable for virtual screening settings with a common binding pocket. We bench- mark our scoring functions on two simplified docking-related tasks: decoy pose scoring and rigid conformer docking. Our method attains similar but faster perfor- mance on crystal structures compared to the widely- used Vina and Gnina scoring functions, and is more robust computationally predicted structures. Code is available at https: / /gggtb...b.ttt//gg///ggg..... com hammmjjingg20 6/scalar- fields.
INTRODUCTION
Proteins are the macromolecular machines that drive almost all biological processes, and much of early-stage drug discovery focuses on finding molecules which bind to and modulate their activity. Molecular docking- the computational task of predicting the binding pose of small molecule to protein target- is important step in this pipeline. Traditionally, molecular docking has been formulated as an optimization problem over a scoring function designed to computational proxy for the free energy (Torres al., 2019; Fan et al., 2019). Such scoring functions are typically a sum of pairwise interaction terms between atoms with physically-inspired functional forms and empirically tuned weights (Quiroga & Villarreal, 2016). While these terms simple and hence fast to evaluate, exhaustive sampling or optimization over the space of ligand poses is difficult and leads to the significant runtime of docking software.
ML-based scoring functions for docking have been an active area of research, ranging in sophistica- tion from random forests to deep neural networks (Yang et al., 2022; Crampon al., 2022). These efforts have largely sought to more accurately model the free energy based on docked pose, which important for downstream identification of binders versus non-binders (virtual screening). How- ever, they have not addressed nor reduced the computational cost required to produce these poses in the first place. Hence, independently of the accuracy of these workflows, molecular docking for large-scale structure-based virtual screening remains computationally challenging, especially with the growing availability of large coliion- databases such as ZINC (Tingle et al., 2023).
In this work, we explore different paradigm and motivation for machine learning scoring functions, with the specific aim of accelerating scoring and optimization ligand poses for high-throughput molecular docking, To do so, we forego the physics-inspired functional form of commonly used scoring functions, and instead frame the problem as that of learning sealar fields independently as- sociated with the 3D structure of the protein and ligand, respectively. We then define the score to be the cross-correlation between the overlapping scalar fields when oriented according to the ligand pose. While seemingly more complex than existing scoring functions, these cross-correlations can

---

Published as a conference paper at ICLR 2024
B
& Φ
ESF
ESF
Figure 1: Overview of the scalar field-based scoring function and docking procedure. The transla- tional FFT procedure is shown here; the rotational FFT is similar, albeit harder to visualize. (A) The protein pocket and ligand conformer are independently passed through equivariant scalar field net- works (ESFs) to produce scalar fields. (B) The fields are cross-correlated to produce heatmaps over ligand translations.
(C) The ligand coordinates are translated to the argmax of the heatmap. Addi- tional scalar field visualizations are in Appendix C.
be rapidly evaluated over a large number of ligand poses simultancously using Fast Fourier Trans- forms (FFT) over both the translational space IR3 and the rotational space 50(3). This property allows for significant speedups in the optimization over these degrees of freedom.
Scalar fields representing molecules in 3D space have been previously parameterized as neural fields (Zhong et al., 2019), i.e., neural networks taking coordinates of query point input and producing field values as output. However, our scalar fields must be defined relative to the ligand (or protein) structure-
as graph embedded in 3D space-and be SE(3) equivariant in order to yield an invariant scoring function. To satisfy these requirements, we introduce equivariant scalar field networks (ESFs), which parameterize the scalar fields with message passing E3NNs (Thomas et al., 2018; Geiger & Smidt, 2022). Unlike neural fields, these networks yield a compact represen- tation of the scalar field in a single forward pass and are automatically SE(3) -equiivrrra
Contrasting with existing MI. scoring functions, the computational cost of our method can be amor- tized at several levels of abstraction, significantly accelerating runtimes for optimized workflows. For example, unlike methods that require one neural network forward pass per pose, our network is evaluated once per protein structure or ligand conformer independently. Post-amortization, we attain translational and rotational optimization runtimes as fast 160 μss and 650 μss respectively, with FFTs. Such throughputs, when combined with effective sampling and optimization, could make docking of very large compound libraries feasible with only modest resources.
Empirically, we evaluate our method on two simplified docking-related tasks: (1) decoy pose scoring and (2) rigid conformer docking. On both tasks, scoring function is competitive with buth-but faster than-the scoring functions of Gnina (Ragoza et al., 2017: McNutt et al., 2021) and Vina (Trott & Olson, 2010) on PDBBind crystal structures and is significantly better on ESMFold structures. We then demonstrate the further ad advantages of runtime amortization on the virtual screening- setup of the PDEIOA test set (Tosstorff et al., 2022), where-since there is only one unique protein structure our method obtains 50x speedup in total inference time at no loss of accuracy.
To summarize, our contributions are:
We are the first to propose learning protein-ligand scoring functions based on cross- correlations of scalar fields in order to accelerate pose optimization in molecular docking.
We formulate a neural network parameterization and training procedure for learning equiv- ariant scalar fields for molecules.
We demonstrate that the performance and runtime of our scoring function holds promise when evaluated on decking-related tasks.
2

---

Published as a conference paper at ICLR 2024
2 BACKGROUND
Molecular docking. The two key components of a molecular docking algorithm are (1) one or more scoring functions for ligand poses, and (2) a search, sampling, or optimization procedure. There is considerable variation in the design these components and how they interact with each other, ranging from exhaustive enumeration and filtering (Shoichet et al., 1992; Meng et al., 1992) to genetic, gradient- based, or MCMC optimization algorithms (Trott& & Olson, 2010; Morris et al., 1998; McNutt al., 2021). We refer to reviews elsewhere (Ferreira al., 2015; Torres al., 2019; Fan et al., 2019) for comprehensive details. These algorithms have undergone decades of develop- ment and have been given rise to eelllestablished software packages in academia and industry, such as AutoDock (Morris & Lim- Wiiny, 2008), Vina (Trott & Olson, 2010) and Glide (Halgren et al., 2004). In many of these, the scoring function is designed not only to identify the binding pose, but also to predict the binding affinity or activity of the ligand (Su et al., 2018). In this work, however, we focus on learning and evaluating scoring functions for the rapid prediction of binding poses.
ML methods in docking. For over decade, ML methods have been extensively explored to im- prove scoring functions for already-docked ligand poses, i.e., for prediction of activity and affinity virtual screens (L.i et al., 2021; Yang et al., 2022; Crampon et al., 2022). On the other hand, developing ML scoring functions as the direct optimization objective has required more care due the enormous number of function evaluations involved. MedusaNet (Jiang et al., 2020) and Gnina (Ragoza et al., 2017; McNutt et al., 2021) proposed to sparsely use CNNs for guidance and re-ranking (respectively) in combination with traditional scoring function. DeepDock (Méndez- Lucio et al., 2021) used hypernetwork to predict complex-specifie parameters of simple statistical potential. Recently, geometric deep learning models have explored entirely different paradigms for docking via direct prediction of the binding pose (Stärk et al., 2022: Zhang et al., 2022; Lu et al., 2022) or via a generative model over ligand poses (Corso et al., 2023).
in structural-basee
FFT methods in docking. Methods based on fast Fourier transforms have been widely applied for the related problem of protein-proter docking. Katchalski-Katte FFTs over the translational space R33 to rapidly evaluate poses using scalar fields that encode the shape complementarity of the two proteins. Later works extended this method to rotational degrees of freedom (Ritchie & Kemp, 2000; Ritchie et al., 2008; Padhorny et al., 2016), electrostatic po- tentials and solvent accessibility (Gabb et al., 1997; Mandell et al., 2001; Chen & Weng, 2002), and data-driven potentials (Neveu et al., 2016). Today, FFT methods are routine step in estab- docking programs such as PIPER (Kozakov et al., 2006), ClusPro (Kozakov et al., 2017), and HDOCK (Yan et al., 2020), where they enable the evaluation of billions of poses, typically as an initial screening step before further evaluation with a more accurate scoring function.
al. (1992) first proposed using
lished protein-prote
In contrast, FFT methods have been significantly less studied for protein-ligand docking. While a few works have explored this direction (Padhorny et al., 2018; Ding et al., 2020; Nguyen et al., 2018), these algorithms have not been widely adopted nor been incorporated into established dock- ing software. A key limitation is that protein-ligand scoring functions are typically more complicated than protein-protein scoring functions and cannot be easily expressed as a cross- cross-correlatin between scalar fields (Ding et al., 2020). To our knowledge, no prior works have explored the possibility of overcoming this limitation by learning cross-correlation based scoring functions.
METHOD
3.1
EQUIVARIANT SCALAR FIELDS
We consider the inputs to a molecular docking problem to be a pair of protein structure and ligand molecule, encoded featurized graphs G G.,, and with the protein structure associated with alpha carbon coordinates X = [xP |xf....x, x Np E 3R3..P. The molecular docking problem is to find the ligand atomic coordinates X [xf,...xx
x NN R33xN, N2 of the true binding pose. To this end, our aim is to parameterize and learn (multi-channel) scalar fields & = 6(x: GP, xP) and &² := (x G , associated with the protein and ligand structures, respectively, such that the scoring function evaluated on any pose X E RR.N, is given by
E(X P x = Σ
/ R3
0e (x GP, x (x; (x,))
(1)
3

---

Published as a conference paper at ICLR 2024
where de refers to the eh channel of the scalar field. While neural fields that directly learn func- tions R + R have been previously developed as encodings of molecular structures (Zhong et al., 2019), such formulation is unsuitable here the field must be defined relative to the variable-sized structure graphs GP, G and transform appropriately with rigid-body motions their coordinates.
Instead, we propose parameterize the scalar field as sum of contributions from each ligand atom or protein alpha-carbon. where each contribution is defined by its coefficients in a spherical harmonic expansion centered that atom (or alpha-carbon) coordinate in 3D space. To do so, we choose a set Ry R+ + IR of radial basis functions (e.g., Gaussian RBFs) in ID and let Y be the real spherical harmonics. Then we define
(x(x;
E n.J.C.m
Aenjtm (G, X)R (Ixx))(((
x Xn X X
where here (and elsewhere) we drop the superscripts L, P for common definitions. Given some constraints how the vector of coefficients Acnjém transforms under SE(3), this parameterization of the scalar field satisfies the following important properties:
Proposition 1. Suppose the scoring function is parameterized as in Equation 2 and for any R E SO(3), t E R3 we have Aonjem(G, R.X. + the (real) Wigner D-matrices, i.e., irreducible representations of SO(3). Then for any g SE),),
= E D Diame (R) (G, (G,.. X) where D (R) are
1. The scalar field transforms equivariantly: (x(x; G. 0g(g-1 -xx
x: G.
2. The scoring function invariant: E(g.XP g..xh) = - (XX.,X).
See Appendix A for the proof. We choose to parameterize A Aenjem (G, R.X) with E3NN graph neural networks (Thomas et al., 2018; Geiger Smidt, 2022), which are specifically designed to satisfy these equivariance properties and produce all coefficients in single forward pass. The core of our method consists of the training of two such equivariant scalar field networks (ESFs), one for the ligand and one for the protein, which then parameterize their respective scalar fields. While the second property (invariance of the scoring function) technically the only one required by the problem symmetries, the first property ensures that different ligand poses related by rigid-body transformations can be evaluated via transformations of the scalar field itself (without re-evaluating the neural network) and is thus essential to method.
Next, we show how this parameterization enables ligand poses related by rigid body motions to some reference pose to be rapidly evaluated with fast Fourier transforms (all derivations in Appendix A). There are actually two ways to do so: we can evaluate the score all poses generated by translations of the reference pose, or via rotations around some fixed point (which always choose to be the center of mass of ligand). These correspond to FFTs over R and SO(3), respectively.
3.2 FFT OVER TRANSLATIONS
We first consider the space of poses generated by translations. Given some reference pose XL, the score as a function of the translation is just the cross-correlation of the fields o² and of::
E _ +t)= xx)pxx++)==>
Σ +
o. (x))))
2x t)
IOo (t)
(3)
where we have dropped the dependence on G, X for cleaner notation and applied Proposition 1. By the convolution theorem, these cross-correlations may be evaluated using Fourier transforms:
1
[
(4)
(2T)3//2
Hence, in order to simultancously evaluate all possible translations of the reference pose, we need to compute the Fourier transforms of the protein and ligand scalar fields. One naive way of doing so would be to explicitly evaluate Equation 2 at an evenly-spaced grid of points spanning the structure and then apply fast Fourier transform. However, this would be too costly, especially during training time. Instead, we observe that the functional form allows us to immediately obtain the Fourier transform via the expansion coefficients Acnjémi
Flo.] (k) Σ c ik-Xn Σ
Σ m,n
Acnjtm (II (k/||k]
(5)
4

---

Published as a conference paper at ICLR 2024
where now y must refer to the complex spherical harmonies and the coefficients must be trans- formed correspondingly, and
2 H H1 V je(kr) dr
the eh order spherical Bessel transform of the radial basis functions. Conceptually, this expression corresponds to first evaluating the Fourier transform of each atom's scalar field in coordinate sys- tem centered at its location, and then translating it via multiplication with e e- in Fourier space. Importantly, He [2y] and Y be precomputed and cached grid of points independently of any specific structure, such that only the translation terms and expansion coefficients need to be computed for every new example.
3.3 FFT OVER ROTATIONS
We next consider the space of poses generated by rotations. Suppose that given some reference pose XL, the protein and ligand scalar fields are both expanded around some common coordinate system origin using the complex spherical harmonics and a set of global radial basis functions Sy (r):
(x(x)=
Σ S(xx (x/[|xl) J.E.m
We seek to simultaneously evaluate the score of poses generated via rigid rotations of the ligand, which (thanks again to Proposition 1) is given by the rotational cross-correlation
E(X", "_,,,XX') =
ESe 6. (x) 'x) d'x -
(8)
Cross-correlations of this form have been previously studied for rapid alignment of crystallographic densities (Kovacs & Wriggers, 2002) and of signals the sphere in astrophysics (Wandelt & Górski, 2001). It turns out that they can also be evaluated in terms of Fourier sums:
J.e de (x) (R xx)
Σ dian Te e
E,m.h.n
where E, 17, w are related to the the Euler angles the rotation R, de is the (constant) Wigner D- matrix for a rotation of n/2 around the y-axis, and
mn Σ BP. J.k
Bejem cken Gjz where Gyk Gz==ʃ S,(r)Sa()) dr
BL
(10)
Thus the main task is to compute the complex coefficients Bejtm of the ligand and protein scalar fields, respectively. This is not immediate as the fields are defined using expansions in "local" radial and spherical harmonic bases, i.e., with respect to the individual atom positions as opposed to the coordinate system origin. Furthermore, since we cannot (in practice) use a complete set of radial or angular basis functions, it is generally not possible to express the ligand or protein scalar field as defined in Equation 2 using the form in Equation 7. Instead, we propose find the coefficients Bejém that give the best approximation to the true scalar fields, in the sense of least squared error.
Specifically, suppose that R E X are the values of Moeal real local basis functions (i.c., different origins, RBFs, and spherical harmonics) evaluated at NNgrd grid points and A E R Noor is the vector of coefficients defining the scalar field de- Similarly define S F NyNss x Nysto using the real versions of the global basis functions. We seek to find the least-squares solution B e R Nawal to the overdetermined system of equations RA SB, which is given by
B = (S"
(11)
Notably, this is simply a linear transformation of the local coefficients Aenjém- Thus, if we can pre- compute the inverse Gram matrix of the global bases (S" -1 and the inner product of the global and local bases S R. then for any new scalar field de the real global coefficients are immediately available via a linear transformation. The desired complex coefficients can then be easily obtained via a change of bases. At first glance, this still appears challenging due the continuous space of possible atomic or alpha-carbon positions, but an appropriate discretization makes the precomputa- tion relatively inexpensive without significant loss of fidelity.
5

---

Published as a conference paper at ICLR 2024
3.4 TRAINING AND INFERENCE
We now study how the rapid cross-correlation procedures presented thus far are used in training and inference. For a given training example with protein structure X the scoring function E(X x²) should ideally attain a maximum at the true ligand pose X X=*. We equate this task to that of learning an energy based model to maximize the log-likelihodd of the true pose under the model likelihood p p(x*) x cxp [E(X* P. x.)]. However, as is typically the case for energy- -aased models, directly optimizing this objective is difficult due to the intractable partition function.
Instead, following Corso et al. (2023), we conceptually decompose the ligand pose X to tuple x4=(x= R. t) consisting of a zero-mean conformer X rotation R. and a translation t. from which the pose coordinates are obtained: X R.X +t. Then consider the following conditional logllikelihoods:
log X2, R) = FX x²) log /
exp [F(XP R.X +t'] d't!
(12a)
R3
log p(R xx,t) = E("", x²) log
exp [x' -t, -R,, dR'
(12b)
50(3)
We observe that these integrands are precisely the cross-correlations in Equations 3 and 8, respec- tively, and can be quickly evaluated and summed for all values of t' and H' using fast Fourier transforms. Thus, the integrals- are the marginal likelihoods p(x², Rt) and p(X², -t-are tractable and the conditional log-likeiihoods can be directly optimized in order to train the neural network. Although neither technically corresponds to the joint lig-liieldhood of the pose, find that these training objectives work well in practice and optimize their sum in our training procedure.
At inference time, a rigid protein structure X is given and the high-level task is score or optimize candidate ligand poses XL_ A large variety of possible workflows can be imagined; however, for proof of concept and for our experiments in Section 4 we describe and focus on the following relatively simple inference workflows (presented in greater detail in Appendix B):
Translational FFT (TF). Given conformer X we conduct grid-based search over R and use FFT to optimize t in order to find the best pose (xC, RR t). To do so, we compute the Fourier coefficients (Equation 5) of the protein X once and for each possible ligand orientation R. We then use translational cross-correlation (Equation 3) to find the best translation t for each R and return the highest scoring combination.
Rotational FFT (RF). Given a conformer X we conduct a grid-based search over t and use FFT to optimize R. To do so, we compute the global expansion coefficients B Bejen of the protein X relative its (zero) center of mass (Equation 11). We then rotational correlaorlssss (Equation 8) to find the best orientation It for each t and return the highest scoring combi- nation.
relative to each possible ligand position t and once for the ligand XC
Translational scoring (TS). Here we instead are given list of poses (X R, t) and wish to score them. Because the values of R nor may not satisfy a grid structure, we cannot use the FFT methods. Nevertheless, we can compute the (translational) Fourier coeffi- cients of the protein X and for each unique oriented conformer R.X of the ligand using Equation 5. We then evaluate
E(XP, R.X +t) -
[lo."(k) Flo ( (.XX))(k)))) e- ik-t d'k
(13)
Since the Fourier transform is an orthogonal operator on functional space, this is equal to the real-space cross-correlation.
score a list of poses (xC, R, t) using the global spherical expansions Bejem- We obtain the real expansion coefficients of the protein relative to each t and for each ligand conformer X using Equation 11. The score for (XC, Rt t) is then given by the rotational cross-correlation
Rotational scoring (RS). Analogously,
E(XP, R. +t) = Σ Berm (XP -t) -t) (XC) (R)G
(14)
C.j.k.E.m...
where Gjk is as defined in Equation 10 and De are the real Wigner D-matrices.

---

Published as a conference paper at ICLR 2024
Table 1: Typical runtimes of the computations involved in inference- time scoring and optimization procedures, measured on PDBBind with one V V100 GPU. The three sets of rows delineate computa- tions that are protein-specifer gggandsspecifin or involve both protein and ligand, respectively.
Frequency
Computation
TF
RF
TS
RS
Runtime
Coefficients Aenjem
v
65 ms
Per protein structure
FFT coefficients
/
V
7.0 ms
+ Per translation
Global expansion Bejém
v
80 ms
Coefficients AenjEm Global expansion Bestm FFT coefficients
V
4.3 ms
Per ligand conformer
17 ms
+ Per rotation
1.6 ms
Per conformer x rotation Translational FFT Rotational FFT Translational scoring Rotational scoring
160
Per conformer x translation
v
650
1.0 js
Per pose
8.2 ps
The runtime of these workflows can vary significantly depending on the parameters, i.e., number of proteins, ligands, conformers, rotations, and translations, with amortizations possible several levels. Table provides a summary of the computations in each workflow, their frequencies, and typical runtimes. We highlight that the RF workflow is well-suited for virtual screening since the precomputations for the protein and ligand translations within a pocket can be amortized across all ligands. Furthermore, if the ligands are drawn from a shared library, their coefficients can also be precomputed independent of any protein, leaving only the rotational FFT as the cost per ligand- protein pair. See Appendix F for further discussion on runtime amortization.
4 EXPERIMENTS
We train and test our model on the PDBBind dataset (Liu et al., 2017) with splits as defined by Stärk al. (2022). We train two variants of our model: ESF and ESF-N, where the latter is trained with rotational and translational noise injected into the examples to increase model robustness. All other model and inference- time hyperparameters are discussed in Appendix E. For the test set, we consider both the co-crystal structures in the PDBBind test split and their counterpart ESMFold complexes as prepared by Corso et (2023). We also collect a test set of 77 crystal structures (none of which are in PDBBind) of phosphodiesterase 10A (PDE10A) with different ligands bound to the same pocket (Tosstorff et al., 2022). This industrially--suuced dataset is representative a real-world use case for molecular docking and benchmarks the benefits of runtime amortization with our approach.
To evaluate our method against baselines, we note that a scoring function by itself is not directly comparable complete docking programs, which also include tightly integrated conformer search, pose clustering, and local refinement algorithms. Here, however, we focus on the development of the scoring function itself independently of these other components. Thus, consider two simplified settings for evaluating our model: (1) scoring decoy poses with the aim of identifying the best pose among them, and (2) docking rigid conformers to given pocket, similar the re-docking setup in Stärk et al. (2022). The first setting focuses on evaluating only the quality of the scoring function itself, whereas the second a simplified version of typical docking setting that circumvents some of the confounding factors while still allowing the benchmarking of FFT-accelerated optimization.
For the baselines, we select Gnina (MeNutt et al., 2021) as the representative traditional docking software, which runs parallel MCMC chains to collect pose candidates that are then refined and re-ranked to produce the final prediction. Among the scoring functions supported by Gnina, we evaluate its namesake CNN (Ragoza et al., 2017) scoring function and the traditional scoring func- tion of Vina (Trott& & Olson, 2010). For recent deep learning baselines, we adapt TANK (Lu et al., 2022) and DiffDook (Corso et al., 2023) for decoy-pose scoring. Specifically, with DiffDook we use the provided confidence podell-wiiee is natively trained as an <2 Á RMSD classifier to score all poses. With TANKBind, we use the predicted distance map and ground truth conformer interatomic distances to calculate Egeneration as the pose score. We do not evaluate these methods on pocket-level conformer docking they cannot be easily adapted for this task.
7

---

Published as a conference paper at ICLR 2024
Table Decoy scoring results (median over test complexes). All RMSDs are heavy-atom symme- try awarc. The best results from our method (ESF) are underlined if not bolded. 'These methods have been adapted for the pose scoring task; see main text.
Crystal structures
ESMFold structures
Time per
<2 Á
Top RMSD
Top
%
Á AUROC
Top RMSD
Top Rank
%
Method
Pose
Complex
AUROC
Rank
<2 Å
<2
Vina
0.93
0.54
2
91
0.86
2.43
419
43
3.4 ms
110 426
Gnina
0.90
0.59
3
83
0.84
2.19
1110
46
13.0 ms
TANKBind
0.69
4.01
6811
10
0.64
4.22
8538
9
1.1
62 ms
DiflDock!
0.96
0.66
3
87
0.89
2.01
143
50
ms
2041 s
ESE-TS
0.87
0.59
3
87
0.82
1.38
24
57
1.0 jus
3.2
ESF-RS
0.87
0.63
3
85
0.82
1.75
22
53
8.2 jus
5.7 8
ESF-N-TS
0.92
0.69
4
81
0.87
1.64
22
54
1.0 js
3.2
ESF-N-RS
0.92
0.75
5
80
0.87
1.74
26
53
8.2 us
5.7
4.1
SCORING DECOYS
For each PDBBind test complex, we generate 323 - = 32767 decoy poses by sampling 31 trans- lational, rotational, and torsional perturbations to the ground truth pose and considering all their possible combinations. On median, the RMSD of the closest decoy is 0.4 Å, and 1.6% of all poses (n 526.5) are below 2 Á RMSD (Appendix D.1). We then score all poses using the baselines and with our method in both TS (Equation 13) and RS (Equation 14) modes. The quality of each scoring function is evaluated with the AUROC when used as <2 RMSD classifier, the RMSD of the top-ranked pose (Top RMSD), the rank of the lowest-RMSD pose (Top Rank), and the fraction of complexes for which the top-ranked pose is under 2 Á RMSD.
As shown in Table 2, our method is competitive with the Gnina and Vina baselines on crystal structures and better on ESMFold structures. This improved robustness is expected since the in- teraction terms in traditional scoring functions are primarily mediated by sidechain atoms, which are imperfectly predicted by ESMFold, whereas our scalar fields only indirectly depend on the sidechains via residue-level coefficients. Our method matches and outperforms the DifiDock confi- dence model on crystal and ESMFold structures, respectively, in terms of identifying the best poses, although DiffDock obtains better <2 À RMSD AUROC. Finally, our method substantially ourper- forms TANKBind in all metrics. Overall, the performance of our method is superior in the TS mode relative RS, likely due the spatially coarser representation of the scalar fields in the global spherical harmonic expansion (i.e., Equation 7) relative to the grid-based Cartesian expansion.
In terms of runtime per pose, our method faster than Vina by several orders of magnitude, with even greater acceleration compared to the CNN-based Gnina and the E3NN- DiffDock, both of which require one neural network evaluation for every pose. The runtime improvement per complex is more tempered since the different proteins and ligand in every complex limit the opportunity for amortization. In fact, of the total runtime per complex in Table 2, only 1% (TS) to 5% (RS) is due to the pose scoring itself, with the rest due to preprocessing that must be done for every and ligand independently. Similar to our method, TANKBind requires only a single neural network evaluation per complex and thus has comparable per-pose runtime, but enjoys faster per-complex runtime due to the lack of such preprocessing.
protein
4.2 DOCKING CONFORMERS
We consider the task of pocket-level docking where all methods are given as input the ground-truth conformer in a random orientation. Following common practice (McNutt et al., 2021), we aim to provide 4 Å translational uncertainty around the true ligand pose in order to define the binding pocket. To do so, we provide Gnina with bounding box with 4 Á of padding around the true pose, and provide our method with a cube of side length Á as the search space for (with a random grid offset). For PDE1OA, and cross-dock to that protein structure. In all docking runs, we deactivate all torsion angles so that Gnina docks the provided conformer to the pocket.
define the pocket using the pose of the first listed complex (PDB 5SFS)
8

---

Published as a conference paper at ICLR 2024
Table 3: Rigid conformer docking results (median over test complexes). All RMSDs are heavy- atom symmetry aware. The median RMSD of our method (ESF) is lower-bounded at 0.5-0.6 A by the resolution of the search grid (Appendix G.2). The runtime shown as an average per complex, excluding/ / including pre-computations that can be amortized.
PDBBind test
Crystal
ESMFold
PDEIOA
%
Med. RMSD
Med. RMSD
%
Med.
Method
<2 Á
<2 A
Runtime
<2 A
RMSD
Runtime
Vina
79
0.32
24
6.1
20
74
0.75
6.1
Gnina
77
0.33
28
5.9
23
73
0.77
6.0
ESF-TF
70
1.13
31
4.6
s/ s/3
67
1.20
1.0
ESF-RF
71
0.97
32
4.4
67
73
0.82
0.5 1.5
ESF-N-TF
72
1.10
2.9
8/
64
1.11
1.0
/
ESF-N-RE 73
1.00
47
3.0
0.5
s/
68
70
1.00
0.5 1.5
As shown in Table 3, the baseline scoring functions obtain excellent success rates on the PDBBind crystal structures (79%). Our method is slightly weaker but also obtains high success rates (73%). The performance decrease in terms of Median RMMS is somewhat larger, likely due to the coarse search grid over non-FFT degrees of freedom (Appendix G.2) and the lack of any refinement steps (which are an integral part of Gnina) in our pipeline. On ESMFold structures, however, our method obtains nearly twice the success rate (47% 28%) of the baseline scoring functions. Unlike in decoy scoring, noisy training noticeably contributes to the performance on ESMFolld structures, and the procedure generally outperforms TF, likely due to the relatively finer effective search grid in rotational cross-correlations (Appendix G.2).
Because of the nature of the PDBBind workflow, the total runtime comparable to or slower than the baselines when precomputations are taken into account. However, in terms of the pose optimization itself, our method significantly faster than the Gnina baselines, despite performing a brute force search over the non-FFT degrees of freedom. While it is also possible to trade-off performance and runtime by changing various Gnina settings from their default values, our method expands the Pareto front currently available with the Gnina pipeline (Appendix G.3; Figure 8). This favorable tradeoff affirms the practical value-add of our method in the context of existing approaches.
To more concretely demonstrate the runtime improvemets: of method with amortization, then dock the conformers the PDE10A dataset. Our method again has similar accuracy to the baselines (Table 3); however, because of the common pocket, all protein-level quantities are com- puted only once and the total runtime is significantly accelerated. For the RF procedure in particular, the computation of protein global coefficients on the translational grid is by far the most expensive step (Table 9), and the remaining ligand precomputations are very cheap. The amortization of these coefficients leads to a 45x speedup the overall runtime (67 s + 1.5 s). (The runtime for Gnina is also accelerated, although to a lesser extent, due to the smaller ligand size.) As the number of ligands increases further, the total runtime per complex of our method would further decrease.
CONCLUSION
We have proposed a machine-learned based scoring function for accelerating pose optimization in molecular docking. Different from existing scoring functions, the score is defined cross- correlation between scalar fields, which enables the use FFTs for rapid search and optimiza- tion. We have formulated novel parameterization for such scalar fields with equivariant neural networks, as well as training and inference procedures with opportunities for significant runtime amortization. Our scoring function shows comparable performance but improved runtime on two simplified docking-related tasks relative to standard optimization procedures and scoring functions. Thus, our methodology holds promise when integrated with other components into full docking pipeline. These integrations may include functions, and architectural adaptations for conformational (i.e., torsional) degrees of freedom all potential directions of future work.
search, refinement with traditional scoring

---

Published as a conference paper at ICLR 2024
ACKNOWLEDGMENTS
We thank Hannes Stärk, Samuel Sledzieski, Ruochi Zhang, Michael Brocidiacono, Gabriele Corso, Xiang Fu, Felix Faltings, Ameya Daigavane, and Mario Geiger for helpful feedback and discussions. This work was supported by the NIH NIGMS under grant # #IR35GM14186) and Department of Energy Computational Science Graduate Fellowship.
REFERENCES
Rong Chen and Zhiping Weng. Docking unbound proteins using shape complementarity, desolva- tion, and electrostatics. Proteins: Structure, Function, and Bioinformatics, 47(3))281-292) 2002.
Gabriele Corso, Hannes Stärk, Bowen Jing, Regina Barzilay, and Tommi Jaakkola. Diffdock: Diffusion steps, twists, and turns for molecular docking. In The Eleventh International Conference on Learning Representations, 2023.
Kevin Crampon, Alexis Giorkallos, Myrtille Deldossi, Stéphanie Baud, and Luiz Angelo Steffenel. Machine-learning methods for ligand proteenn molecular docking. Drug discovery today, 27(1): 151-164, 2022.
Xinqiang Ding, Yujin Wu, Yanming Wang, Jonah Z Vilseck, and Charles L Brooks III. Accelerated cdocker with gpus, parallel simulated annealing, and fast fourier transforms. Journal of chemical theory and computation, 16(6):3910-3919 2020.
Jiyu Fan, Ailing Fu, and Le Zhang. Progress in molecular docking. Quantitative Biology, 7:83-89, 2019.
Leonardo G Ferreira, Ricardo N Dos Santos, Glaucius Oliva, and Adriano 1) Andricopulo. Molecu- lar docking and structure-band based drug design strategies. Molecules, 20(7):13384-13222 2015.
Henry A Gabb, Richard M Jackson, and Michael JE Stembeerg. Modelling protein docking using shape complementarity electrostatics and biochemical information. Journal of molecular biology, 172(1):106-120, 1997.
Mario Geiger and Tess Smidt. e3nn: Euclidean neural networks. arXiv preprint arXiv:2207.09453, 2022.
Thomas A Halgren, Robert B Murphy, Richard A Friesner, Hege S Beard, Leah L Frye, W Thomas Pollard, and Jay L Banks. Glide: a new approach for rapid, accurate docking and scoring. 2. enrichment factors in database screening. Journal of medicinal chemistry, 2004.
Huaipan Jiang, Mengran Fan, Jian Wang, Anup Sarma, Shruti Mohanty, Nikolay V Dokholyan, Mehrdad Mahdavi, and Mahmut T Kandemir. Guiding conventional protein-ligand docking soft- ware with convolutional neural networks. Journal of Chemical Information and Modeling, 60 (10):4594-4602, 2020.
Bowen Jing, Gabriele Corso, Jeffrey Chang, Regina Barzilay, and Tommi Jaakkola. Torsional diffu- sion for molecular conformer generation. arXiv preprint arXiv:2206.01729, 2022.
Ephraim Katchalski-Katzin Isaac Shariv, Miriam Eisenstein, Asher A Friesem, Claude Aflalo, and Ilya A Vakser. Molecular surface recognition: determination of geometric fit between proteins and their ligands by correlation techniques. Proceedings of the National Academy of Sciences, 89 (6):2195 2199, 1992.
Julio A Kovacs and Willy Wriggers. Fast rotational matching. Acta Crystallographica Section D: Biological Crystallography, 1888 1286, 2002.
Dima Kozakov, Ryan Brenke, Stephen R Comeau, and Sandor Vajda. Piper: docking program with pairwise potentials. Proteins: Structure, Function, and Bioinformatics, 65 (2):392-000 2006.
fft-based protein
10

---

Published as a conference paper at ICLR 2024
Dima Kozakov, David R Hall, Bing Xia, Kathryn A Porter, Dzmitry Padhorny, Christine Yueh, Dmitri Beglov, and Sandor Vajda. The cluspro web server for protein-protere docking. Nature protocols, 12(2):255-22 20088,8,8888 2017.
Hongjian Li, Kam-Heung Sze, Gang Lu, and Pedro J Ballester. Machine-learning scoring functions for structure-based virtual screening. Wiley Interdisciplinary Reviews: Computational Molecular Science, 11(1)ce1478 2021.
Zhihai Liu, Minyi Su, Li Han, Jie Liu, Qifan Yang, Yan Li, and Renxiao Wang. Forging the basis for developing protein-ligand interaction scoring functions. Accounts of Chemical Research, 50 (2):302 10002-30,,000 2017.
Wei Lu, Qifeng Wu, Jixian Zhang, Jiahua Rao, Chengtao Li, and Shuangjia Zheng. Tankbind: Trigonometry-awaa neural networks for drug-protens binding structure prediction. Advances in neural information processing systems, 2022.
Jeffrey G Mandell, Victoria A Roberts, Michael E Pique, Vladimir Kotlovyi, Julie Mitchell, Erik Nelson, Igor Tsigelny, and Lynn F Ten Eyck. Protein docking using continuum electrostatics and geometric fit. Protein engineering, 14(2):005-113, 2001.
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, 2021.
Oscar Méndez-Lucio, Mazen Ahmad, Ehecatl Antonio del Rio-Chanona, and Jörg Kurt Wegner. A geometric deep learning approach predict binding conformations of bioactive molecules. Nature Machine Intelligence, 3112):0033-1339 2021.
Elaine C Meng, Brian K Shoichet, and Irwin D Kuntz. Automated docking with grid-based energy evaluation. Journal of computational chemistry, 13(4):505-52) 5244, 1992.
Garrett M Morris and Marguerita Lim- Lim-Wilby. Molecular docking. In Molecular modeling of pro-
teins, pp. 365-382 Springer, 2008.
Garrett M Morris, David S Goodsell, Robert S Halliday, Ruth Huey, William E Hart, Richard K Belew, and Arthur J Olson. Automated docking using a lamarckian genetic algorithm and an empirical binding free energy function. Journal of computational chemistry, 19(14):169--166, 1998.
Emilic Neveu, David W Ritchie, Petr Popov, and Sergei Grudinin. Pepsi-dock: a detailed data-driven protein interaction potential accelerated by polar fourier correlation. Bioinformatics, 32 (17):1693-i701, 2016.
Trung Hai Nguyen, Huan-Xiang Zhou, and David DL Minh. Using the fast fourier transform in binding free energy calculations. Journal of computational chemistry, 39(1 1):6216636, 2018.
Dzmitry Padhorny, Andrey Kazennov, Brandon S Zerbe, Kathryn Porter, Bing Xia, Scott Mottarella, Yaroslav Kholodov, David W Ritchie, Sandor Vajda, and Dima Kozakov. Protein- protein docking by fast generalized fourier transforms on 5d rotational manifolds. Proceedings of the National Academy of Sciences, 113(30):E4286-BA) E4293, 2016.
Dzmitry Padhorny, David R Hall, Hanieh Mirzaei, Artem B Mamonov, Mohammad Moghadasi, Andrey Alekseenko, Dmitri Beglov, and Dima Kozakov. Protein ligand docking using fft based sampling: D3r case study. Journal of computer-aided molecular design, 32:225-2321 2018.
Rodrigo Quiroga and Marcos A Villarreal. Vinardo: A scoring function based on autodock vina improves scoring, docking, and virtual screening. PloS one, 11(5):c0155183, 2016.
Matthew Ragoza, Joshua Hochuli, Elisa Idrobo, Jocelyn Sunseri, and David Ryan Koes. Protein- ligand scoring with convolutional neural networks. Journal of chemical information and model- ing, 57(4):942-95 9274229957,00 20,,-957,2177, 2017.
David W Ritchie and Graham JL Kemp. Protein docking using spherical polar fourier correlations. Proteins: Structure, Function, and Bioinformatics, 39(2): 2000.
11

---

Published as a conference paper at ICLR 2024
David W Ritchie, Dima Kozakov, and Sandor Vajda. Accelerating and focusing protein preeeni docking correlations using multi dineesinnmnii rotational fft generating functions. Bioinformatics, 24(17):1865-1873, 2008.
Brian K Shoichet, Irwin D Kuntz, and Dale L Bodian. Molecular docking using shape descriptors. Journal of computational chemistry, 13(3):300-377, 1992.
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-20521. PMLR, 2022.
Minyi Su, Qifan Yang, Yu Du, Guoqin Feng. Zhihai Liu, Yan Li, and Renxiao Wang. Comparative assessment of scoring functions: the casf-2016 update. Journal of chemical information and modeling, 59(2):895-99) 993,999933999 2018.
Nathaniel Thomas, Tess Smidt, Steven Kearnes, Lusann Yang, Li Li, Kai Kohlhoff, and Patrick Riley. Tensor field networks: Rotation-and canslation-equivarian neural networks for 3d point clouds. arXiv preprint, 2018.
Benjamin Tingle, Khanh G Tang, Mar Castanon, John Gutierrez, Munkhzul Khurelbaatar, Chin- zorig Dandarchuluun, Yurii Moroz, and John J Irwin. Zinc-22- free database of tangible compounds for ligand discovery. Journal of Chemical Information and Mod- eling, 63(4): 166-11666661 1176, 2023.
Pedro HM Torres, Ana CR Sodero, Paula Jofily, and Floriano P Silva-Jr. Key topics in molecular docking for drug design. International journal of molecular sciences, 20(18):4574, 2019.
Andreas Tosstorff, Markus G Rudolph, Jason C Cole, Michael Reutlinger, Christian Kramer, Hervé Schaffhauser, Agnès Nilly, Alexander Flohr, and Bernd Kuhn. high quality, industrial data set for binding affinity prediction: performance comparison in different early drug discovery scenarios. Journal of Computer-Aided Molecular Design, 36(10):753-76,, 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):55--461 2010.
Benjamin D Wandelt and Krzysztof M Górski. Fast convolution on the sphere. Physical review D, 63(12):123002, 2001.
Wikipedia. Hankel transform Wikipedia, the free encyclopedia, 2023. URL https:
//en. wwkipedia htapppeiaa.rr//eeeeeeeeeee ciieesioimessions.
Yumeng Yan, Huanyu Tao, Jiahua He, and Sheng- You Huang. The hdock server for integrated protein-protein docking. Nature protocols, 15(5):182--1852. 2020.
Chao Yang, Eric Anthony Chen, and Yingkai Zhang. Protein-ligant docking in the machine- learning era. Molecules, 27(14):4568, 2022.
Anna Yershova, Swati Jain, Steven M Lavalle, and Julie C Mitchell. Generating uniform incremental grids on so (3) using the hopf fibration. The International journal of robotics research, 29(7):801- 812, 2010.
Yangtian Zhang, Huiyu Cai, Chence Shi, Bozitao Zhong, and Jian Tang. E3bind: An end-to-end equivariant network for protein-ligand docking. arXiv preprint 1t.v:22100060100 2022.
Ellen D Zhong. Tristan Bepler, Joseph H Davis, and Bonnie Berger. Reconstructing continuous dis- tributions of 3d protein structure from cryo-em images. In International Conference on Learning Representations, 2019.
12

---

Published as a conference paper at ICLR 2024
A MATHEMATICAL DETAILS
A. PROOF OF PROPOSITION 1
Proposition Suppose the scoring function parameterized as in Equation 2 and for any R S0(3), E R we have A R.X. +t) = Emn the (real) Wigner D-matrices, L.e,, irreducible representations of SO(3). Then for any g E SE(3),
mm (R) Menjem((G,)) X) where (R(R) are
1. The scalar field transforms equivariantly: de 6x g.).X) 6c 1 xx; G.
2. The scoring function is invariant: EXg.XP, 9 F(XP, X4).
Proof. Let the action of g (R, E SE(3) be written Rd [x - t). We first note that X g.x. 1 X Xn || and R8 (x g.xn) = 9 x Xn. Then
+ Rxx + t and hence g x +
Ox(x; G,9.X)= = E
Aonjem (G, R.X t) (1xx g.xx
x g-Xnn
|| - g.x
n.j.l.m
Σ
Aenzem (G,X)R, g.x..ll Σ D (R)Y (
x 9.Xnn
||x g.Xn
m.j.d.mm
RF"(x - g.8n) Aenjém ( (Gx I.x (
Σ
||x 11x-9-0 9.X
n.j.tim
9 g-1,x x-11 Xnn
= Σ
Aengem (G x - xx, 11
(1g-1 - X,
n.j.d.mm
6elg Xx;G G)
Next,
Elgxx ,gxx")
&e
(x
g. (x; G¹ ,g. dx
R3
Σ ( x: G. 20 x; at xt) dxx
Σ/ (x': x 1ot x d'x
Σ
where the last line has substitution x' =g'xx 'x with g volume preserving on R3
A.2 DERIVATIONS
In this section we describe the derivations for the various equations presented in the main text. We use the following convention for the (one-dimensiona)) Fourier transform and its inverse:
Flf((k) = 22T f f-tha ike f(x) dx
(15a)
[ (ikzz [(k) dk
(15b)
2n
Equation It well known (Wikipedia, 2023) that given function over R with complex spher- ical harmonic expansion
f(r)
Σ Sem (11
(16)
t.m
its Fourier transform is given by
f(k) Σ (-1) (
(17)
E,m
13

---

Published as a conference paper at ICLR 2024
where
Fr,m (k) V. L.
dr
(18)
with Jr the -order Bessel function of the first kind. Relating these to the spherical Bessel functions je via J+41/2(z) = 22x/ mrj(z)) we obtain
/ π
2 /-- fe.m fe.m(r)je(kr)+ r² dr mm
Fm (k)
(19)
which is the form of Equation 6. To apply this to our scalar fields, we define the translation operator Tx [xx = f fx r) and note its composition with the Fourier transform
(Fo Tf] -rr Fifl
(20)
We then decompose the form of our scalar fields (Equation 2) into contributions from zero-origin spherical harmonic expansions
6x(x)
Σ Tx., [den](x))
(21a)
den (x)
Σ f,m
Σ
Rensem Yi" (x/[|>]])
(21b)
3
Penti
Hence, the Fourier transform of each contribution is
Flden](k/||
(--)) ((kk))) (k/[[]]
(22)
f,m
Equation 5 is then obtained via Equation 20 and the linearity of the Fourier and spherical Bessel transforms.
Equation 9 We source (with some modifications) the derivation from Kovacs & Wriggers (2002). We consider the cross-correlation
of dR) = / d(x)v(R ddx
(23)
which is the same as Equation with o of and v = of since of is real field. Expanding in complex spherical harmonics Y and radial bases Sy:
((x)
Σ
(xx) = Σ V Vean (IxI)
(24)
j.
j.cm
We then obtain
c(R)
Σ
De Pame R))2yen pemm
[S, -X YX
dx
(25a)
d.P mmm.
Σ
Dc aa (R) j8m po
[S) Syl(r) (2 dr
[Ym
Y (f) df (25b)
5.3 L. mannn
52
Gzy
dzc 8mn
Σ (R)
D mm
Σ j.j
4 jem ye G
(25c)
Now to evaluate the complex Wigner D-matrix, we adopt the extrinsic zyz convention Euler angles (applied right-to-left) and note that any rotation (0,0, 2b) be decomposed
R( 0, v = R2 - (2)) (T R
7/2)
(26)
Next, can easily check (using the standard spherical harmonics) that the Wigner D-matrix for rotation about the z-axis is diagonal and given by Dr (Rz (R)(w))= = 8mne- innw Hence,
D mn R( f x)) = e -ime d - -hu "h e -kon
(27)
where de (R(RR((T(2)) constant and real. Complex conjugation then gives Equation
14

---

Published as a conference paper at ICLR 2024
Equation 12
The conditional likelihood is
logp(t x R) log p(x², ,,,))
(28a)
p(xc, .R)
= logp(XX .R,t) - log
p(x R.,t)
(28b)
Rd
= log E(X", x! - log
exp
[E(X", R.X +t)] d'
(28c)
R3
Similarly,
logp(I xtt))= log p(X², t,t)
(29a)
p(Xc,t)
= loogp(X. BR,t) log
p(x ,,,t)
(29b)
50(3)
log FX( x4) log
exp [F(XP R.X +t)] dR'
(29c)
50(3)
Finally, we move to the protein coordinates (invoking the invariance of the score E) to obtain form consistent with the rotational cross-correlations (Equation 8).
Equation 13 Given pose X = =R.X. + t, we evaluate
E(X", X.X +t)= /. (xx +t)dxx
(30)
The functional inner product is equivalent in Fourier space:
F(X", R.X +t) +t)= Σ/ Flo?((k) (: + t)](k) d
(31)
Then with the translation operator T defined previously,
(x; R.X +t) = To(o( Rxxx)
(32a)
FloE(:: R.X + t))(k) = e-kk ( ()()
(32b)
We then substitute into Equation 31 to obtain Equation 13.
Equation 14 Given pose X R.X + we assume that the field (-;'" t) and o x) are written in the real global spherical harmonic expansion:
&e X
- t) = - Ber S3 (](1 (Xx])
(33a)
j.f.m
05 xx: x
Σ
Begem By (x]
(33b)
j.l.
Then, analogously to Equation 25,
F(X", C.X +t) F(X" R,,X
(34a)
= =/ R3
X"
-t)o (R- Ix; xO) d'x
(34b)
Drmm (R)
Σ
Blsam Bly Gjy
(34c)
m.
341"
Complex conjugation has been omitted because the coefficients and D-functions are now real.
15

---

Published as a conference paper at ICLR 2024
B ALGORITHMIC DETAILS
Below, we present detail the four inference procedures introduced in Section 3.4. The three blocks of computations are color-coded corresponding to protein preprocessing (green), ligand pre- processing (blue), and the core computation (red) and labelled with typical runtimes from Table (unlabelled lines have negligible runtime). The various loop levels make clear that depending on the workflow, the protein and ligand processing precomputations can be amortized and approaches a negligible fraction of the total runtime. Note, however, that for readability we have presented the algorithms assuming that all possible combinations (i.e., of proteins, ligand conformers, rota- tions, and translations) are of interest; if this not true (for example in PDBBind, in any typical pose-scoring setting), then the full benefits of amortization may not be fully realized.
Algorithm 1: TRANSLATIONAL FFT
Input: Proteins ( conformers {(G), x
Output: Docked poses (X, X Vi
foreach (G.². X do
// protein preprocessing
Compute coefficients A = {Ap cjnem (G.P x.)} with neural network; // 65 ms // 7.0
Compute Fourier-space field values F( using x
foreach (G x do
// ligand preprocessing
Compute coefficients {A , } with neural network; // 4.3 ms
foreach R3 [gid C S3(3) do
Compute rotated coefficients A A, using D' (RR))
Compute Fourier-space field values Fi using A k X
// 1.6
foreach
do foreach (. do foreach R2 E [R] gnid C S3(3) do 1,,Ke, VT using FFT;
// pose zation
Compute F(X," P. + +t,
11 160 μ
Ei th + {max, arg max E( R + t)
k* + arg max, E;: X - Rs x + ti.:;
Algorithm 2: ROTATIONAL FFT
Input: Proteins {(G],, X,"), }} conformers ((G), X/1)
Output: Docked poses (X." X Vi
foreach (G. X do
// protein preprocessing
Compute coefficients A, {AP ( (Gf X.)) with neural network; // 65
for ty (t) C 3 do
L Compute global expansion B.e (Bettm) from Af. X to // 80
foreach (Gh, x do Compute coefficients A {A (G), x()) with neural network; // 4.3 ms
// ligand preprocessing
Compute global expansion B/d (Bosmm) from A x
// 17 ms
foreach (G, x.") do foreach (G x do
1/ pose opt imi maization
foreach t, [t] C R do
Compute F.X.
tr. R.X(),, YR using FFT; Ei, R2 +- {max, arg max) EX. tR, R.X/ 4 + t)
1/ 650 μt
ks + arg max Et: X/A t R / tg+
16

---

Published as a conference paper at ICLR 2024
Algorithm 3: TRANSLATIONAL SCORING
Input:
Proteins {(G} )) conformers {(G x rotations {Re}, translations (tz) Output: Scores E R + tt) Vi,h, k,
foreach (G ( do
// protein preprocessing
Compute coefficients A, (A cjnem (G.. XF)} with neural network; // 65 ms
Compute Fourier-space field values Fi using A x
// 7.0
foreach (G/ X do
// ligand preprocessing
Compute coefficients A = [AA (G } with neural network;
// 4.3 ms
foreach R2 do
Compute rotated coefficients R using D' (R²))
Compute Fourier-space field values ]i using A X
// 1.6
foreach
do foreach (Gt, x do
// scoring
foreach R2 do
foreach t
Compute E(X!", R2 X + t²) using Equation 13:
// 1.0 μs
Algorithm 4: ROTATIONAL SCORING
Input: Proteins {(G), X()),, conformers {(G)) X rotations (Ra}, translations (tz) Output: Scores EXX," R, + +te) Vi, k,
foreach (G. ( do
// protein preprocessing Compute coefficients A, = {AP (Agnee (G.. x()) with neural network; // 65 ms
for ty E {t) Jgid C 3 do
L Compute global expansion B. [Bejtm] from A, x, - th // 80 ms
foreach (G/,, x do Compute coefficients A/ {AL (jem (G/ ( with neural network; // 4.3 // 17 ms
// ligand preprocessing
Compute global expansion B/C (Bojmm) from A x
foreach (G.², X.")
// scoring
foreach (G x do
foreach R2 do
foreach t; do Compute E(X,, R2 X +t) using Equation 14;
// 8.2 μs5
17

---

Published as a conference paper at ICLR 2024
C
LEARNED SCALAR FIELDS
6qqw
6jap
6np2
Figure 2: Visualizations of learned scalar fields. All five channels of the ESF-N learned scalar fields d² (top row) and 6º (bottom row) are shown on the ry-plane passing through the center of mass of the ligand, with box diameter of 20 A. Positive values of the field are in blue and negative values in red. At left, the ligand and pocket structures are shown looking down the z-axis. Note that as the fields are only 2D slices, not all 3D features visible the structures are visible in the fields.
18

---

Published as a conference paper at ICLR 2024
6uvp
dua
6oxq
6jsn
6moa
Figure 3: Visualizations of learned scalar fields, continued.
19

---

Published as a conference paper at ICLR 2024
D EXPERIMENTAL DETAILS
D.1
DECOY SET
Given a zero-mean ground-truth ligand pose X we generate 323 32767 decoy poses via the following procedure.
Sample 31 translational pertubations: (V, = .1..31 and set to 0, with units in A.
Sample 31 rotational perturbations: Rj FromRotvec(ts). 0.500. 1...31 and set R I3.
Sample 31 noisy conformers X = .1. by sampling torsional updates Atk (0, (T/2)Lmm where Nr is wrapped normal distribution (Jing et al., 2022) and m the number of torsion angles. The torsional updates are applied to the smaller side of the molecule. Set X = -X"
Set X ijk = R t, j,j,j 0...3I and discard X X
PDB ID 6A73 is excluded from the procedure due to the high level of graph symmetry and significant runtime for computing RMSDs for all decoys. Summary statistics for the decoy sets of the remaining 362 PDB IDs are presented in Figure 4.
1.0
30000
25000
20000
# of decoys
CDF
0
5
10
15
20
15000
0.4
0.002
50
10000
0.2
5000
0.000
0.0
0
0.0
0.2
0.4
0.6
0.8
1.0
5
10
15
20
RMSD
Figure 4: Decoy set statistics. Top left: histogram of RMSDs across all decoys sets (12 M total). Botton left: histogram of minimum RMSDs among the decoy sets. All sets have pose less than RMSD À from the true pose. Right: cumulative density function of RMSDs each decoy set. Bottom right inset: all decoy sets have at least 23 poses with RMSD <2 A.
D.2 RUNTIME MEASUREMENTS
All runtime measurements were performed on a machine with Intel Xeon Gold 6130 CPUs and 8 Nvidia Tesla V100 GPUs. Gnina was run with default thread count settings. All of our processes were run on single V100 GPU. For our method, we performed runtime analysis using CUDA events to remove the effects of asynchronous CUDA execution. Script loading, model loading, and algorithmic-level precomputations (which, if necessary, can cached on disk) were excluded from the analysis. For Gnina, we attempted to remove similar overhead by timing single-poss scoring- only runs as representative of constant overhead costs. We report conformer docking runtimes in Table using the PDBBind crystal structure;; ESMFold runtimes marginally shorter. Typical
20

---

Published as a conference paper at ICLR 2024
runtimes reported in Table 1 and Appendix B are obtained from timing runs with our method across the entire PDBBind test set.
D.3 DATASETS
As noted previously, we use train, validation, and test splits from Stärk et al. (2022). However, due to RDKit parsing issues with Gnina-docked poses, the following 30 complexes are excluded (leaving 333 remaining) from all rigid conformer docking comparisons against Gnina: 6HZB, 6E4C, 6PKA, 6E3P, 6OXT, 6OY0, 6HZA, 6E6W, 6OXX, 6HZD, 6K05, 6NRH, 60XW, 6RTN, 6D37, 6HLE, 6PY0, 60XS, 6E30, 6HZC, 6Q38, 6E7M, 60IE, 6D3Y, 6D40, 6UHU, 6CJP, 6E3N, 6Q4Q, 6D3X. Scoring comparisons include all test complexes except 6A73, for which decoy poses could not be generated.
We download the 77 PDB IDs provided in Tosstorff et al. (2022) from the PDB to form the PDEIOA dataset, keeping the A chain of each assymetric unit and the Ligand of Interest (LOI) interacting with it. We then align all ligands the crystal structure of 5SFS using the procedure described in Corso et al. (2023) for aligning ESMFold structures, except transforming the ligand rather than the protein. This constitutes the construction of cross-docking dataset due to the of the same pocket for all ligands. Due to RDKit parsing errors with the Gnina-docked poses, the following PDB IDs are excluded from all comparisons: 5SFA, 5SED, 5SFO, 5SEV, 5SF9, 5SDX, SSFC. The remaining 70 ligands are shown superimposed on the 5SFS pocket in Figure 5.
Figure 5: PDE10A ligands aligned on 5SFS
21

---

Published as a conference paper at ICLR 2024
E HYPERPARAMETERS
Model hyperparameters The main model hyperparameters are the aspects of the scalar field pa- rameterization (Equation which dictate the number of coefficients Aenjém predicted by the neural network. Intuitively, these determine the expressivity or "intrinsic resolution" of the scalar field which can be parameterized by the ESF neural network. These hyperparameters and their search spaces are detailed in Table 4, with full supporting experimental results in Appendix G.I.
Our data featurization and E3NN neural network architecture are largely adapted from the default settings of Corso et al. (2023). Both the protein and ligand neural networks operate on all heavy atom nodes; however, in the protein network only the alpha-carbons nodes emit Acnjén coefficients that contribute to the scalar field. For simplicity, we also omit ESM features.
Table 4: Training- saainigg-time model hyperparameters and their search spaces. The parameters correspond- ing to experiments the main text are in bold.
Parameter
Search Space
# scalar field channels radial basis functions
3, 5, 8
RBF cutoff radius
3 A, 5 A, Á
Order of spherical harmonics
1,
Inference hyperparameters The primary inference-time hyperparameters control the "search res- olution" over the FFT and grid-search degrees of freedom. For the grid-search, the resolution is con- trolled directly by the selection of grid points. For the FFT, the resolution is effectively controlled by frequency domain representation of the scalar field. Specifically,
In the translational case, this resolution is expressed terms a sampling resolution which corresponds the grid of frequency vectors at which the translational Fourier co- efficients (Equation 5) are evaluated. To generate these frequency vectors, we use a 40 Á cubical domain with periodic boundary conditions.
In the rotational case, this resolution is given by the maximum order of spherical har- monics in the global spherical basis (Equation 7) used for evaluating rotational cross- correlations.
The search space for these resolution hyperpameters is detailed in Table 5. Increasing the search resolution improves accuracy but at the cost of core optimization runtime. The default settings chosen for the main experiments are based qualitative assessment of the best tradeoff. Figure visualizes this tradeoff, with full supporting experimental results Appendix G.2.
Table Inference-time model hyperparameters and their search spaces. The search grid points over SO(3) generated according to Zhong et al. (2019); Yershova al. (2010). The translational search grid points are defined by evenly spacing the stated number of points along each dimension of 8 A cube. The parameters corresponding to experiments in the main text are in bold.
Procedure
Parameter
Search Space
Translational signal resolution Rotataional search grid points
0.5 A, 1 576, 4608
TF
Rotational signal order Translational search grid points
10, 25, 50 73, 93, 133
RF
22

---

Published as a conference paper at ICLR 2024
Translational Optimization (TF)
Ligand preprocessing runtime (ms)
Cross-correlation runtime (ms)
Median RMSD (A)
1.1
7196
715
grid pts.
576
1.53
1.5
930
928
100
122
0.5
i
0.5
Translational cross-correlation resolution (A)
Rotational Optimization (RF)
Protein preprocessing runtime (seconds)
Pross-correlation runtime (ms)
Median RMSD (A)
13a
1.05
0.9
197
200
994
1430
Trans. grid pts.
ga
1.16
1.0
0.98
64
63
332
476
1486
73
1.25
1.15
1.16
30
30
31
157
224
704
Rotational cross-correlation resolution (fmmx)
Figure 6: Inference-time hyperparameters and their effect on performance and average runtime. In the TF procedure (top), the protein preprocessing time is roughly constant and is not shown; similarly in the RF procedure (bottom) the ligand preprcessing time constant. The selected default settings are boxed. All numbers taken from full results PDBBind in Appendix G.2
23

---

Published as a conference paper at ICLR 2024
F RUNTIME AMORTIZATION
In this section we discuss in greater detail the practical implications of runtime amortization for various common workflows molecular docking, focusing on rigid conformer docking as proxy lask. Conceptually, the total runtime of a docking workflow is given by
( unique) (proteins)
(# unique ligand)
Opeot
+ +Cig
+ Cpaic (# complexes)
(35)
conformers
where Cpror. Clig, Caar refer to the protein precomputations, ligand precomputations, and cross- correlation (i.c., pose optimization) runtimes, respectively. If the number of complexes grows faster than either the number of unique proteins or ligand conformers individually, then the runtime con- tributions from the latter be amortized across the complexes-meaning that the unit cost per complex decreases as their number increases.
To illustrate this amortization, suppose that we have N protein pockets (for simplicity we assume each is on a seperate protein) and M rigid conformers and wish to dock all MN complexes. Then, following the exposition and typical runtimes provided in Section 3.4, compute (Table 6) the total docking cost of this workflow for the TF procedure:
(72 ms) M + (7400 ms) N + (740 ms) MN
(36)
Similarly, the total cost under the RF procedure is:
(58400 ms)- M + (21 ms) N (470 ms) MN
(37)
Note that Gpon < Gig for the TF procedure, whereas Cpre >> Cit for the RF procedure. This is be- cause that the grid search over rotations in TF involves a new set of translational Fourier coefficients for each rotation of the ligand, while the protein is kept fixed; over translations in RF requires a new global spherical expansion Bajam of the protein around each grid point, but only once for the ligand around its center of mass. This discrepancy has significant implications for the suitability each procedure depending on the behavior M, M:
the other hand, the grid search
In the virtual screening- like setting, we keep a fixed protein (M 1) and let the number of ligands increase (N - oo). Hence, the TF procedure is unsuitable, as the cost per ligand (and thus per complex) is constant z seconds, approximately on the same order as traditional docking methods. On the other hand, the cost per ligand under the RF procedure begins at nearly z=60 s with a single ligand, but asymptotically approaches R550 ms as the number of ligands increases (Figure 7, left).
In the inverse-screening like setting, we keep fixed ligand (N = 1) and let the number of proteins increase (M >> oo). Now the RF procedure is unsuitable, as the cost per protein is a constant ==00 which is far too high be useful. On the under hand, the cost per protein under the TF procedure begins at &88 and asymptotically approaches x800 ms as the number of proteins increases (Figure 7, right).
In the general setting, we divide both Equation 36 and 37 by MN to find that when M >> (i.e., more proteins than ligands), the TF procedure is more efficient, and when M & N (more ligands than proteins), the RF procedure is more efficient. However, both procedures converge to similar per-complex runtimes (740 ms vs 470 ms) in the limit of
M, +
24

---

Published as a conference paper at ICLR 2024
Table 6: Computing the runtime of the TF and RF procedures as function of the number of proteins and ligands N. and are the number of points used the grid search over rotations and translations for TF and RF, respectively; they fixed constants the order of ==100 (we use = 729, R = 4608 by default). The rows in red contribute the greatest preprocessing runtime.
Translational Optimization (TF)
Rotational Optimization (RF)
Computation
Count Runtime
Computation
Count Runtime
Protein scalar field
M x 65 ms
Protein scalar field
Mx
65 ms
Protein FFT coeffs
M 7.0 ms
Protein global coeffs
T TM
80 ms
Ligand scalar field
Nx 4.3 ms
Ligand scalar field
N X
4.3 ms
Ligand FFT coeff's
RNx
1.6 ms
Ligand global coeffs
Nx
17 ms
Cross-correlation
RMN 160 js
Cross-correlation
TMNX
650 ps
(Equation 36)
(Equation 37)
Virtual screening (M 1)
Inverse screening (W=1) 1)
60
60
Runtime per complex (s)
TF RF
50
50
40
30
30
20
10
10
100
10
102
103
100
101
102
N (# ligands)
M (# proteins)
Figure 7: Amortized conformer docking runtimes per protein-ligand complex in the virtual screening (left) and inverse screening (right) settings.
comple
25

---

Published as a conference paper at ICLR 2024
G ADDITIONAL RESULTS
G.1 MODEL. HYPERPARAMETER
In Table 7, we explore the impact varying model hyperparameters that determine the expressivity of the scalar field from their default values. We use the ESF-N model and evaluate on the rigid conformer docking task on PDBBind crystal structures. As the results illustrate, increasing the expressivity of the scalar field can bring additional improvements in performance. On the other hand, decreasing the expressivity degrades the performance our method to extent that is no longer competitive with the baselines (Table 3).
Table 7: Rigid conformer docking results for varying model hyperparameters. Runtimes are shown averaged over complexes, excluding including pre pre-computations that can be amortized.
TF
RF
%
Med. RMSD
Med. RMSD
Method
Runtime
Runtime
<2 Á
<2 Á
Baseline
72
1.10
s
73
1.00
0.5 68
# channels = 3
70
1.07
71
0.98
0.4 70
# channels
71
1.15
0.9
1.09
0.7 s/
# RBFs 3
63
1.27
0
64
1.13
0.5 s
RBFs =
77
1.07
0.3s
0.95
s/
RBF cutoff A
72
1.16
s
1.00
s/56/68
RBF cutoff Å
76
1.16
/
71
1.08
0/5
Spherical harmonics max =
62
1.35
0.96
62
1.23
0.5 / 62
Spherical harmonics max
80
0.99
0 10s
0.87
0.5 s/
G.2 INFERENCE HYPERPARAMETERS
In Tables 8-11 below, explore the impact of inference-time hyperparameters on the performance and runtime of our method on the rigid conformer docking task. We use the ESF-N model variant and experiment with the PDBBind crystal test set and PDE10A test set. As discussed in Appendix E, for the TF procedure, we adjust (1) the number of grid points over 50(3) (576 or 4608); and (2) the spatial resolution for translational cross- cross-correlation (either adjust (1) the number of translational grid points in the 8 Á cube (73, 93, or 139); and (2) the resolution for the rotational cross-correlation, expressed in terms of the maximal spherical harmonic order (10, 25, 50). The rows corresponding results in the main Table 3 are bolded.
or 0.5 Å). For the RF procedure,
In all rows, the effective number of poses searched over via both degrees of freedom is computed. To provide an idea of the impact of discretization, we compute the median RMSD of the closest grid point to the ground-truth pose (decomposed into rotational and translational contributions). This serves as a hard lower bound for the median RMSD the output docked pose. In the TF procedure, increasing the resolution is memory intensivee: thus, the RF procedure is more effective leveraging FFT to conduct fine-grad fine-grained search over the accelerated degree freedom. The default reported performance is attained with a translational offset of 0..4 and a rotational offset of 0.16 A. While performance improves with smaller grid offsets, the returns are rapidly dimishing.
The runtime of the method (averaged over 333 PDBBind complexes and 70 PDEIOA complexes) reported and color-coded according to Appendix B: protein preprocessing (green), ligand prepro- cessing (blue), and the pose optimization (red). The effect of the non-FFT grid resolution is also color-coded, i.c., in TF the explicit enumeration over SO(3) grid points directly scales the ligand preprocessing, whereas in RF the cnumeration over H scales the protein preprocessing. As the tables show, the preprocessing of these explicit grid points contributes to the majority of the non- amortizeable runtime. In general, the S0(3) grid ligand preprocessing in TF is less expensive, however, it cannot be amortized when moving from PDBBind to PDEIOA (where the ligands are still distinct). On the other hand, the g grid protein preprocessing time in RF is significantly reduced (very roughly on the order of 70-fold, as expected) in PDEIOA compared to PDBBind.
26

---

Published as a conference paper at ICLR 2024
Table 8: PDBBind TF
Grid offset
Runtime (ms)
Trans. resol.
SO(3) grid pts.
Effective poses
Med. RMSD
%
Prot.
Lig. prep.
Tr.
Rot.
All
Opt.
<2
prep.
1 Á
576
420k
0.52
0.84
0.98
1.53
63
65
931
100
1 1
4608
3.4M
0.50
0.42
0.67
1.10
72
72
7196
715
0.5 Á
576
2.8M
0.25
0.80
0.84
1.50
64
70
928
123
Table 9: PDBBind RF
Grid offset
Runtime (ms)
Trans. grid pts.
5O(3)
Effective poses
Med. RMSD
%
Prot.
Lig. prep.
Tr.
Rot.
All
Opt.
mmax
<2
prep.
73
10
3.2M
0.65
0.38
0.80
1.25
70
30k
85
158
73
25
45M
0.67
0.15
0.70
1.15
69
31k
87
225
73
353M
0.65
0.08
0.67
1.16
70
32k
704
g3
10
6.8M
0.49
0.36
0.64
1.16
73
64k
85
333
9a
25
97M
0.50
0.15
0.53
1.00
73
67k
87
476
93
751M
0.51
0.08
0.52
0.98
74
63k
84
1487
133
10
20M
0.33
0.37
0.51
1.05
74
198k
85
995
133
25
291M
0.33
0.15
0.37
0.90
72
200k
86
1430
Table 10: PDEIOA TF
Grid offset
Runtime (ms)
Trans. resol.
S3) grid pts.
Effective poses
Med. RMSD
%
Prot.
Lig. prep.
Tr.
Rot.
All
Opt.
<2
prep.
I A
576
420k
0.51
0.88
1.00
1.85
56
22
761
89
1 A
4608
3.4M
0.50
0.48
0.69
1.11
64
21
6159
736
0.5 Á
576
2.8M
0.26
0.89
0.93
2.05
50
20
756
106
0.5 Á
4608
23M
0.26
0.44
0.51
1.00
73
20
6147
892
Table 11: PDE10A RF
Grid offset
Runtime (ms)
Trans. grid pts.
SO(3)
Effective poses
Med. RMSD
%
Prot.
Lig. prep.
Tr.
Rot.
All
Opt.
mmax
<2
prep.
73
10
3.2M
0.72
0.38
0.83
1.60
54
476
44
161
73
25
45M
0.57
0.16
0.59
1.21
63
549
42
227
73
50
353M
0.65
0.08
0.65
1.30
64
635
59
718
93
10
6.8M
0.46
0.39
0.63
1.05
64
1014
42
327
93
25
97M
0.48
0.16
0.51
1.00
70
946
43
465
g3
751M
0.49
0.09
0.50
0.99
64
943
1483
133
10
20M
0.34
0.41
0.55
1.17
64
2798
42
986
133
25
291M
0.33
0.16
0.36
0.96
69
2912
45
1469
27

---

Published as a conference paper at ICLR 2024
G.3 FERFORMANCE-RUNTE TRADEOFF
In Figure 8, we further investigate the tradeoff between speed and performance offered by our method compared to Gnina (with the Vina scoring function). While in the main results (Ta- ble 3) we run Gnina using all default settings, it is possible to reduce the runtime (and perfor- mance) by adjusting these settings. In particular, we explore setting --max _onaa_m_______ _cc_seps and --minimizee_iie iters to 5 independently and in combination. Together with the default runs and the --ssorre_nnl runs, these trace out a Pareto frontier representing the tradeoff between runtime per complex and <2 Å RMSD success rate. With the default settings, Gnina outperforms all variants of our method on the PDBBind crystal and PDEIOA test sets. However, Figure shows that we can reach previously inaccessible regions in the accuracy v.s. runtime tradeoff landscape.
0.8
0.7
RMSD
0.6
0.5
% <
A
0.4
0.4
0.3
0.5
0.2
0.2
0.0
0.1
4
0.0
0.0
5
10
15
20
4
PDBBind Runtime (s)
PDE10A Runtime (s)
Figure 8: Tradeoff between speed and accuracy using our method compared to Gnina on PDBBind crystal structures (left) and PDEIOA (right). In both cases, variants of our method (blue dots) enable possibilities not reachable with Gnina (orange curve).
28


Reviews:
[{"review_id": "NGQbyr0EOF", "reviewer": "Reviewer_ncny", "summary": "This paper explores the use of machine learning to accelerate the process of molecular docking, specifically by learning a scoring function that allows for more rapid optimization over rigid-body degrees of freedom using fast Fourier transforms. The scoring function is defined as the cross-correlation of multi-channel ligand and protein scalar fields, and is parameterized by equivariant graph neural networks. The method is benchmarked on two simplified docking-related tasks: decoy pose scoring and rigid conformer docking. The authors demonstrate that their method achieves similar performance to existing scoring functions but with faster runtime. They also highlight the advantages of runtime amortization in virtual screening settings with a common binding pocket.", "questions": "- The evaluation of the proposed method is limited to simplified docking-related tasks, such as decoy pose scoring and rigid conformer docking. While the results on these tasks are promising, it is unclear how the method would perform in more complex and realistic docking scenarios.\n\n- The paper compares against Vina and Gnina scoring functions in experiments. I would like to know how this method compares with other SOTA methods/scoring functions [1, 2, 3, 4].\n\n[1] Corso, Gabriele, et al. \"Diffdock: Diffusion steps, twists, and turns for molecular docking.\" arXiv preprint arXiv:2210.01776 (2022).\n\n[2] Stärk, Hannes, et al. \"Equibind: Geometric deep learning for drug binding structure prediction.\" International conference on machine learning. PMLR, 2022.\n\n[3] Koes, David Ryan, Matthew P. Baumgartner, and Carlos J. Camacho. \"Lessons learned in empirical scoring with smina from the CSAR 2011 benchmarking exercise.\" Journal of chemical information and modeling 53.8 (2013): 1893-1904.\n\n[4] Lu, Wei, et al. \"Tankbind: Trigonometry-aware neural networks for drug-protein binding structure prediction.\" Advances in neural information processing systems 35 (2022): 7236-7249.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "- The paper proposes a novel approach for learning protein-ligand scoring functions based on cross-correlations of scalar fields. This approach enables the use of Fast Fourier Transforms (FFTs) for rapid search and optimization, which can significantly accelerate pose optimization in molecular docking.\n\n- The paper introduces runtime amortization techniques that lead to significant speed improvements in the scoring function. By precomputing certain coefficients and leveraging the common pocket structure in docking, the total runtime of the method is accelerated, making it more efficient for large-scale docking tasks.\n\n- The proposed scoring function shows comparable performance to existing scoring functions on two simplified docking-related tasks: decoy pose scoring and rigid conformer docking. In some cases, the proposed method even outperforms existing scoring functions, particularly on ESMFold structures.", "weaknesses": "- Lack of Detailed Hyperparameter Analysis: The paper mentions that hyperparameters are detailed in Appendix D, but the specific hyperparameters used in the experiments are not provided in the extracted content. Without this information, it is difficult to understand the impact of different hyperparameter settings on the performance of the proposed method.\n\n- Limited Discussion on Runtime Amortization: The paper mentions the benefits of runtime amortization with the proposed method, but there is limited discussion on the practical implications and potential limitations of this approach. Further analysis and discussion on the tradeoff between runtime improvements and accuracy would provide a more comprehensive understanding of the method's performance.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "SdDBQ9HbRw", "reviewer": "Reviewer_BpA5", "summary": "This paper proposes learning scalar fields for proteins and ligands using equivariant neural networks. The fields are cross-correlated to define a scoring function that enables rapid optimization over rigid-body degrees of freedom with fast Fourier transforms (FFTs).", "questions": "1. Are there any ablation experiments on FFT？\n2. Why is there no significant improvement in the scoring of crystal structures?", "limitations": "", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. The idea of learning cross-correlation based scoring functions to enable FFT optimization is interesting.\n2. The experiments demonstrate the speed advantage over traditional tools", "weaknesses": "1. Too many model settings to choose from in practice.\n2. The method is only evaluated on simplified tasks, not a complete docking benchmark. Performance in a full virtual screening workflow remains to be tested.\n3. Lack of baselines.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "yUPtlJ2pfI", "reviewer": "Reviewer_9kCY", "summary": "The paper discusses a machine learning approach to speed up molecular docking in virtual screening. The authors introduce a scoring function using fast Fourier transforms and cross-correlation of ligand and protein scalar fields, proving more efficient on benchmark tasks than prior methods. It performs quicker on crystal structures and more reliably on computationally predicted structures.", "questions": "- Could you please clarify what type of information is captured by the different channels?\n- Is the method predominantly dependent on geometric data, such as shape complementarity? Is the chemical data thoroughly utilized?\n- Ellen (cryoDRGN2) lays out a coarse-to-fine approach that subdivides the grid into finer segments. Would it be feasible to incorporate such a technique into this project?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "- The speedup mainly comes from factorization and amortization. Factorizing the interaction into two scalar fields is a novel and effective idea. By computing cross-correlation, the cost of pre-computations can be amortized over a group of operators.\n- The manuscript is well organized, clearly written and easy to follow.\n- The experimental results are impressive.", "weaknesses": "- The selection of hyper-parameters seems opaque in this document. There's no clear indication of most hyper-parameters being fine-tuned, or it's not showcased, leaving readers unsure about which hyper-parameter has the most substantial influence on the final outcome. For instance, how were decisions made regarding the number of channels or the quantity of spherical harmonics?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)