For NeuralDock inputs, we used crystal structures of 3875 known protein-ligand pairs from the PDBbind 2017 refined set (Wang et al., 2004), which were shuffled into a training set (N = 2,712, 70% of structures), and a validation set (N = 1,163, 30% of structures). One of these structures was not processable by MedusaDock 2.0, and 127 structures exceeded our computational resources. A further 496 ligands were not processable by rdkit 2020.09.3²⁷. This left a training set of N = 2,279 (84.0% of original) and a validation set of N = 972 (83.6% of original). The core set of PDBbind 2013 was used as the test set (N = 195); 154 of these were processable by rdkit. There is no overlap of proteins among the training, validation, and test sets.

We extracted the atoms in the protein which were within a cube of side length 20 angstroms, centered at the ligand. These protein pockets were encoded as 10 × 10 × 10, 2-angstrom resolution images with 8 channels corresponding to a one-hot encoding of (no atom, C, O, N, S, P, H, and other). If multiple atoms were contained in the same 2-angstrom cube, we took the maximum of each channel of the one-hot encodings. Using rdkit, we encoded the ligand as a length 36 atom type vectors with 7 channels (no atom, C, N, O, F, S, and other) and a 36 × 36 bond adjacency matrix with 5 channels (no bond, single, double, triple, aromatic/conjugated). For ligands with greater than 36 heavy atoms (hydrogens were excluded), we removed atoms with the least bond order until we were left with 36 total heavy atoms, thus attempting to preserve the important topologies present in the ligand. These input dimensions were chosen so that a comparable number of parameters in NeuralDock would be devoted to processing the protein and the ligand each (Figure 1A), as well as to prevent the ligand representation from becoming too sparse. We also chose a limit of 36 heavy atoms since molecules with greater than 36 heavy atoms are likely to exceed the 500 Dalton cutoff for Lipinski’s rule of five.

For each protein-ligand pair, we ran MedusaDock for 24 h or 1000 iterations (whichever came first) on a single core of an Intel Xeon E5-2,680 v3 processor with 6 GB RAM. We collected summary statistics (mean, median, standard deviation, minimum, maximum, skew, and kurtosis) on the 13 interaction energies computed in MedusaDock’s force field, MedusaScore (Yin et al., 2008). Five of these energies were zero for most or all structures (see Supplementary Information). An interaction energy is defined as the total energy of the protein-ligand complex E P − L minus the contributions from the protein E P and the ligand E L when they are isolated. E int = E P − L − E P − E L

In developing our models, we chose to focus on E_{without VDWR}, the interaction energy excluding repulsive van der Waals forces. E_{without VDWR} is the output of MedusaScore which is known to be the most highly correlated with experimental binding affinity (Yin et al., 2008).

The compounds (N = 997,402,117) were downloaded from the ZINC database (Irwin and Shoichet, 2005) (available in Tranches) and processed into the tensor input format described above. Of these compounds, 936,054,166 (94%) were processable by rdkit. We chose the 1UXM SOD1 PDB structure as a protein target (Berman et al., 2000). The interface between the A and B chains (Figure 2) was chosen as the binding pocket for 1UXM.

For NeuralDock, we chose a fully connected architecture, with spectral normalization (Miyato et al., 2018), dropout (Srivastava et al., 2014) with rate 0.2, LeakyReLU activation (Maas et al., 2013), and skip connections (He et al., 2016). These characteristics were chosen as standard methods of model regularization to prevent overfitting and vanishing or exploding gradients during training.

The NeuralDock architecture (Figure 1A) was implemented in TensorFlow 2.4.0 (Abadi et al., 2016) and Python 3.7. The inputs were the protein pocket and ligand topology, and the outputs were the 7 summary statistics of the 13 interaction energies computed by MedusaDock as well as the binding affinity pK = log₁₀ K (dissociation/inhibitor constant K_D/I with units of molar). We considered K_D and K_I to be the same for the training purposes. We chose 10 fully connected (FC) hidden layers for each of the three parts (protein encoder, ligand encoder, affinity predictor) of the network (Figure 1A), resulting in approximately 45 million trainable parameters. We varied the number of hidden layers as well as their widths during hyperparameter optimization (Table 1). The loss function was the L² (squared difference) loss between the NeuralDock output energies and the MedusaDock output energies, as well as pKs. The Adam optimizer was used (Kingma and Ba, 2015), with a learning rate of 10⁻⁶. Training of each model took place on one NVIDIA Tesla T4 GPU, and the models were trained to convergence within a week. We trained a convolutional architecture, in which the FC blocks in the protein encoder (top left of Figure 1A) were replaced by spectrally normalized 3D inception modules (Szegedy et al., 2015) with a comparable number of trainable parameters (Table 1).

We used Lipinski’s rule of five (Lipinski et al., 2012) to evaluate small molecule lead quality and drug likeness, as well as the Quantitative Estimate of Drug likeness (QED) (Bickerton et al., 2012). Octanol-water partition coefficients (log P) were extracted from HTML files of the ZINC database, while all other quantities were computed using rdkit (RDKit, 2016).

Least squares regression was performed in Python 3.7 using SciPy 1.6.0 (Virtanen et al., 2020). Analysis of covariance (ANCOVA) was performed in Python 3.7 using Pingouin 0.3.12 (Vallat, 2018).

NeuralDock training on MedusaDock energies and experimental pKs converged, achieving agreement with experimental binding affinities for the test set, the PDBbind 2013 core set (Figure 1B), and the validation set taken from PDBbind 2017 (Figure 1C). As discussed by Francoeur et al. (2020), testing on the core set may offer a biased evaluation of binding affinity prediction performance. We provide the core set correlation for comparison with other methods which solely report that data (Table 2). Given our relatively small training set (2,331 structures) and the massive diversity of potential protein pockets and ligands, NeuralDock was able to learn the binding affinities accurately. We believe that the key to NeuralDock’s success is using high quality data produced by MedusaDock, which sampled thousands of conformations for each protein-ligand pair, while supplying coarse 3D protein information and only the small molecule’s topology. By using only the topology of the small molecule as input, we forced NeuralDock to approximate the effects of conformational sampling.

To test the robustness, generalizability, and speed of NeuralDock, we performed virtual screening of a massive library of ligands (N = 936,054,166) against the pocket at the dimeric interface of SOD1 (Figure 2). We compared relationships among MedusaDock E_{without VDWR}, NeuralDock E_{without VDWR}, and experimental binding affinities (pK), and found that they agree (Figure 3). E_{without VDWR} is the component of MedusaScore which is most highly correlated with experimental pK (Yin et al., 2008). We demonstrate the correlation between NeuralDock and MedusaDock E_{without VDWR}, first on the validation set drawn from PDBbind 2017 refined set (r = 0.83, p < 0.0001), and then on 100 random small molecules from ZINC docked to 1UXM, an A4V mutant of the human superoxide dismutase-1 enzyme (r = 0.69, p < 0.0001) (Figure 3A). Therefore, NeuralDock was successful in learning to predict the minimum E_{without VDWR} in a single shot. Additionally, we performed 2-way ANCOVA, which measures the effects of a categorical variable. In this instance, it showed no statistically significant difference in the trends for the validation set and the external validation on 1UXM (F = 0.67, p = 0.41), which provides evidence that NeuralDock successfully generalizes from the training set to other protein-small molecule pairs, as well as achieving cross-docking success. Since the native structure 1UXM was not bound to any small molecule and MedusaDock performs fully flexible conformational sampling of both the protein and the small molecule, our results provide external validation of NeuralDock and demonstrate robustness under cross-docking, i.e., docking of small molecules to a protein as it appears in its native, unbound conformation.

We found agreement of MedusaDock E_{without VDWR} and NeuralDock E_{without VDWR} with experimental pK, with similar correlations for both tools (r = −0.48 for both, p < 0.0001), and no statistically significant difference in correlations with 2-way ANCOVA (F = 1.27, p = 0.26) (Figure 3B). The correlation of MedusaDock E_{without VDWR} with experimental pK (r = −0.48) is comparable to that of AutoDock Vina scoring (r = 0.41) as reported in Francoeur et al. (2020), with MedusaDock performing better, likely due to our extensive sampling and computational effort for each protein-ligand pair. NeuralDock predicts experimental pK better than MedusaDock on the validation set (Figures 1B, C, 3B), however, experimental validation of the pKs is needed to confirm that this result holds for SOD1 and other proteins. NeuralDock’s agreement with MedusaDock predictions demonstrates that it may be useful in replacing traditional docking tools such as MedusaDock and AutoDock Vina.

For the billion-molecule docking of 1UXM, we find that the 100 small molecules with the highest NeuralDock-predicted pKs also have low binding energies, with agreement of NeuralDock and MedusaDock E_{without VDWR} on those small molecules (Figure 3C). Furthermore, NeuralDock predicts that few other molecules have such high pK and low energy (Figure 3D). Combining these results with our validation set regression (Figure 1C), we predict that some of the 100 small molecules chosen have pK of approximately 8 (K ≈ 10 nanomolar). Ninety-five of the 100 compounds satisfy Lipinski’s rule of five, and all have reasonable QED scores (median 0.50, range 0.33–0.66). The drug likeness of these compounds can be explained as bias from the ZINC Tranches, in which only 14,744,513 (1.5% of the over 997 million) compounds have molecular weight exceeding 500 Daltons or log P exceeding 5.

Various architectures for NeuralDock yielded similar results (Table 1). We optimized the NeuralDock architecture based on E_{without VDWR} correlation, as E_{without VDWR} is most highly correlated with binding affinity (Yin et al., 2008). Optimizing the network based on E_{without VDWR} prevents overfitting the hyperparameters (i.e., number of layers and total number of trainable parameters), as MedusaDock E_{without VDWR} is an external source of training information. Including other MedusaDock energies in training also helped to prevent overfitting (Supplementary Table S2). Even a limited network such as the 4.9 million parameter model was able to generalize from the training set to the validation set (Table 1).

Benchmarking on a single Tesla T4 GPU was able to predict energies and pKs of 96,000 protein-small molecule pairs in 203.8 s, or 2 ms per protein-small molecule pair. Docking of the 937 million ZINC compounds with NeuralDock took 21 h on 25 GPUs, which matches our benchmarking of 2 ms per protein-small molecule pair per GPU. For 3875 structures (including those that exceeded our computational resources), it took 120 processors 4 weeks to perform extensive MedusaDock sampling (1000 iterations), or an average of 20 h per structure per processor. Note that the ability of MedusaDock to find low energy binding poses requires repeated conformational sampling, which can be controlled by running many iterations of MedusaDock, with each iteration taking seconds to minutes. Taking 10 h as a conservative minimum compute time for 1000 iterations of MedusaDock, NeuralDock performs 1.7 × 10⁷ times faster. NeuralDock performs 10⁵ times faster than DOCK 6⁶, and 100 times faster than K_DEEP (Jiménez et al., 2018). NeuralDock is also faster than the Def2018 General Ensemble by Francoeur et al. (2020), since we directly predict binding affinity in a single forward pass instead of requiring many conformational samples and forward passes through the network (Francoeur et al., 2020).

There are several potential improvements to NeuralDock. First, NeuralDock takes a coarse atomic image of the 20 Å protein pocket as input. Including inputs such as protein amino acid composition, secondary structure, and hydrophobicity may improve NeuralDock predictions. This potential change in inputs is supported by the fact that NeuralDock was able to learn local contributions to the energy well (i.e., E_VWDA, Supplementary Tables S1, S2), while having difficulty with global contributions (e.g., E_solv, Supplementary Tables S1, S2). Second, NeuralDock is trained on a relatively small training set, and we did not use any small molecules which are off-target, i.e., do not bind to the given protein. However, our external validation with MedusaDock indicates that this potential source of bias was at least partially addressed by the variety of binding affinities present in our training set. This deficiency may be remedied by docking more targets with MedusaDock as a reference, such as the PDBbind general set, which incurs a significant computational cost in data preparation. Additionally, the PDBbind general set experimental pKs may be noisier than in the refined set. Third, the binding energies and pK are invariant under rotations of the protein pocket as well as under permutations in the atom order for the small molecule bond adjacency matrix and atom type vector. Newer neural network architectures which directly leverage these physical symmetries have been proposed (Finzi et al., 2020; Wu et al., 2021).

We posit that our small training set is not necessarily a major impediment to NeuralDock performance. The so-called double descent phenomenon occurs when the number of tunable parameters in a statistical model significantly exceeds the number of training samples (or degrees of freedom) (Belkin et al., 2019). The traditional tradeoff between bias and variance constitutes the first descent in test error, in which the optimal number of model parameters lies below the number necessary to precisely fit the data, a number known as the interpolation threshold. Contrary to conventional statistical wisdom, however, over-parametrizing may provide benefit in the model performance, causing a second descent in test error. In certain cases, additional training data may actually hurt model performance, whereas training past the point of 0 training set error (another instance of the interpolation threshold) provides a similar benefit to over-parametrization (Nakkiran et al., 2020). Although the double descent phenomenon has yet to be formally proven for deep neural networks, we observed the phenomenon in practice while training NeuralDock. With 46 million parameters and roughly 2000 data points for training, NeuralDock takes advantage of the second descent to achieve its high accuracy.

MolGAN is a generative adversarial approach to drug discovery, which does not take into account the protein pocket (Goodfellow et al., 2014; de Cao and Kipf, 2018). Since NeuralDock is end-to-end differentiable, one can augment MolGAN training with NeuralDock acting as a scoring function to create drugs targeting specific proteins. The addition of such guidance may improve the stability of MolGAN training. While automatically differentiable force fields are currently being developed for molecular dynamics (Schoenholz and Cubuk, 2019), the advantage of NeuralDock is it can immediately evaluate the quality of a generated small molecule, implicitly performing conformational sampling, whereas conformational sampling must still be done for a differentiable force field. NeuralDock does not predict the explicit minimum energy configuration for the protein-small molecule pair, however, once specific compounds are identified via NeuralDock, the candidates can be docked by MedusaDock or another docking software prior to experimental validation, as the computational cost of doing so is much less than the cost of traditional docking for the entire compound library.