Guiding Conventional Protein–Ligand Docking Software with Convolutional Neural Networks

Data Set.

We train our CNN model based on a refined set of the PDBbind database.⁵⁵ This data set is further curated by Wang et al.¹⁶ The complexes with more than one ligand and the complex with less than two rotatable bonds were removed since those complexes lack flexibility in the ligand. After the compilation, this data set contains 3875 protein–ligand complexes with their experimentally solved 3D structures.

To feed the CNN with data for training, we first obtain the massive center of the binding site by calculating the average coordinates across all atoms of the ligand in the native state complex. After that, we crop a cube of 20 × 20 × 20 ų around the binding site (First step in Figure 1). All of the atoms (both protein atoms and ligand atoms) inside this cube area are retained for training. In the default case, the input grid has the resolution of 20 × 20 × 20, which means each voxel represents the atoms inside that 1 × 1 × 1 ų area.

To train our CNN model, we randomly select 3100 out of the 3875 proteins as training data and the remaining 775 proteins as testing data. For each protein, the data consist of three parts: (1) the ground-truth data, (2) the positive data, and (3) the negative data. The ground-truth data are the X-ray crystallographic structure of the protein–ligand complexes. For positive and negative data, we run MedusaDock for 11 attempts with different random seeds, which generate a total of 230 826 poses for training and a total of 56 098 poses for testing. We consider a pose “positive” if its RMSD is less than 2 Å from the crystallographic structure (ground-truth pose). In our computational experiments, we only include the atom types used by MedusaDock, and we have 10 distinct atom types among all generated poses in the output of MedusaDock. We distinguish each atom type from the protein and the ligand except the hydrogen. The metals are not included in our atom types. We list all atom types in Table S1. As one can expect, each training/testing data is formed as a four-dimensional (4D) tensor, where the first three dimensions indicate the spatial information (X, Y, Z) and the last dimension indicates the atom type in such small cube (1 × 1 × 1 ų) area.

Input Grid Generation.

The 3D grid in Figure 1 is usually implemented by a 3D array, and the value of each element of the array indicates how many atoms are there in the corresponding cube of the 3D grid. As an example, let us assume that we have a 3D grid that represents the space between [−10 Å, −10 Å, −10 Å] and [10 Å, 10 Å, 10 Å] with the origin at the massive center of the binding site. We first translate the 3D space into 3D grid with the resolution of [20 × 20 × 20], where each voxel represents a volume of [1 Å × 1 Å × 1 Å]. Subsequently, we place all of the atoms within such 3D space into the 3D grid based on the 3D coordinate of the atom. More specifically, if we have an atom with the coordinate of (0.2 Å, 3.4 Å, 1.5 Å), we first calculate the index of this atom in the 3D grid (which is [10, 13, 11]). After that, we set the value of the 3D array as 1 at that entry (at the index of [10, 13, 11]). Additionally, using a 4D array, one can distinguish the types of atoms in the protein–ligand complex. The fourth dimension indicates the type of the atoms in that cube.

3D CNN Model.

Traditional energy-based docking approaches suffer from both the accuracy and efficiency perspectives. Although they usually propose several poses for each protein as the candidates based on some scoring functions, however, for different protein–ligand complexes, the threshold of the scoring function for determining whether a pose is good or not is different. Hence, to find the best binding poses, docking attempts need to be performed for a large amount of times or even exhaustively, to search all possible binding poses in the searching area to select promising best poses. It is to be noted that selecting a few poses among thousands of proposed poses is intractable. On the other hand, a CNN model can learn the patterns from the 3D grid of the poses and predict whether each pose is good or bad. Recent studies^47–51 reveal that CNN can play a good role in this binary classification task.

Our CNN model (MedusaNet) in Figure 2 takes the 3D grid of a protein–ligand complex as the input. The first several convolution layers extract the features about the bond and other interactions between nearby atoms. The following fully connected layers reorganize the features and calculate the possibility of whether this input binding pose can be a stable binding. The first convolution layer has 5 × 5 × 5 filters, while all other convolution layers contain 3 × 3 × 3 filters. This is because the first layer requires a large-scale filter to collect more spatial information. The fully connected layers have a size of 1024, and the last layer has two neurons, which collectively capture the possibility of binding or not binding. RELU activation is applied on each convolution layer and fully connected layers. A dropout layer with a dropout ratio of 0.5 is applied after each fully connected layer to reduce the probability of overfitting. The output of the CNN is a value between 0 and 1, which indicates the probability that the pose is a good binding pose. We set the cutoff threshold to 0.5, which means that the poses with output probability less than 0.5 are considered as “bad poses,” whereas others are considered to be “good poses.” Compared to previous neural networks, our proposed CNN model does not employ the pooling layers to avoid the resolution loss during the training.⁵⁷ We also employ MedusaDock energy score as a feature in our CNN model and accept it in the second fully connected layer. This is because the convolution operators cannot catch the features with high-dimensional computation or exponential calculation, while the MedusaDock energy score includes such nonlinear features.

MedusaDock.

MedusaDock starts with a stochastic rotamer library of ligands generation, then shifts and rotates the rotamers into different poses. Based on the Monte Carlo algorithm and score function, MedusaDock tries to search for the best-fit poses within the docking boundary box. During this search process, both protein side chain and ligand will be repacked to imitate the realistic protein–ligand docking procedure. Additionally, to optimize the searching algorithm and reduce the computation of score function, MedusaDock separates the search step into coarse docking and fine docking. In coarse docking, MedusaDock quickly identifies several good poses, and those poses are carefully checked during the fine docking step. To search the docking poses of a protein–ligand complex with MedusaDock, we have to keep running MedusaDock with different random seeds until we cannot find a pose with lower energy. For some proteins, we need to run MedusaDock for thousands of times, even if some good poses have already been found.

Using CNN to Improve MedusaDock.

In general, to solve the protein–ligand docking problems with CNNs, the first step is to convert the protein and ligand molecule file into a 3D coordinate representation. In other words, the 3D coordinate of each atom should be translated to the shift of a 3D grid, where each atom is placed into the corresponding cube. After the 3D/4D tensor is generated as the input for the CNN, the network applies multiple convolution layers to calculate the energy of the bonds and other interactions between the atoms. This information can be captured by the convolution operations and extracted as features. The features are accordingly employed to predict the stability of the pose by the fully connected layers. Recall that, since the 3D/4D array has the position information of each and every atom, the pose of the protein–ligand is already implied from the input tensor.

Although previous CNN-based docking work can achieve successful accuracy in the docking stability prediction problem, to our knowledge, none of the existing approaches can generate the best binding poses directly. To solve the pose prediction problem for a protein–ligand complex, some previous works use traditional energy-based docking software^58–60 to sample a large number of possible binding poses and employ CNN to predict, for each pose, whether it is a stable docking. As a result, such CNN-based approaches cannot improve the docking process explicitly. Hence, in this paper, in addition to proposing a CNN model (MedusaNet), we want to boost the execution and improve the accuracy of the pose prediction with our CNN model.

In the original MedusaDock, to find the best binding pose for a protein–ligand complex, we run MedusaDock for multiple iterations until the global minimum energy is converged.^15–18 More specifically, we record the energy for the output pose from MedusaDock and update the global minimum energy value if we find a pose with lower energy in some iteration. We keep running MedusaDock for more iterations with different random seeds until we cannot find a lower energy pose. It is believed that we have found the best binding pose when the global minimum energy has converged. This might take up to 2000 MedusaDock runs, which corresponds to extremely long execution times, for some complex proteins. To reduce the execution time, we early stop the execution process by checking the output poses with our CNN model. As shown in Figure 2, after each iteration of a MedusaDock run, we send all of the candidate poses to MedusaNet, which in turn predicts whether those poses are good docking or not. Note that the search process early terminates when MedusaNet claims that k good poses are found by MedusaDock, instead of running until global minimum energy converges. We show in the Results section that this early termination strategy reduces the execution time significantly while finding good binding poses for more proteins.

Evaluating the Scoring Capability of the CNN Model.

For this experiment, we compile a data set composed of 3875 protein–ligand complexes from PDBbind v_2016 by retaining the complexes that include only one ligand. Subsequently, the 3875 proteins are randomly divided into training (3100 proteins) and testing (775 proteins) sets. For each protein, we use MedusaDock to perform 11 docking attempts with 11 different random seeds. MedusaDock generates ~286 K candidate binding poses for those proteins. We calculate the root-mean-square deviation (RMSD) of all candidate poses with the native binding pose. If a pose has an RMSD value less than 2 Å, we consider it as a good pose; otherwise, we consider it as a bad pose. We train the CNN model with proteins and ligands in the native pose, good poses, and bad poses in the training set and evaluate the model with the poses of the proteins in the testing set. We applied the cross-validation and randomly split the data set for 10 times to repeat the experiment for 10 times. For all results below, we show the average, minimum, and maximum results across all 10 validations.

We show the performance of our CNN model (MedusaNet) in Table 1 and compare it against three other approaches (MedusaDock energy score, a random forest (RF) approach with 100 decision trees, and another CNN model AtomNet⁴⁸). We compare MedusaNet with a RF approach to check if any simple machine learning method can improve the accuracy. Additionally, we compare MedusaNet against AtomNet to check if our model performs better than other CNN models. AtomNet is an existing CNN model for predicting the binding stability. It takes the 3D grid (with a resolution of 20 × 20 × 20) as the input feature map, and then applies four convolution layers with the parameters of 128 × 5³, 256 × 3³, 256 × 3³, 256 × 3³, respectively (number of filters × filter size). After that, the feature map is sent to two fully connected layers with the size of 1024. We use three different evaluation metrics—accuracy, AUC, and average RMSD—in our experiments. Accuracy is calculated by dividing the total number of poses by the sum of true positives and true negatives. It indicates how many poses we predict correctly. AUC is the area under receiver-operating characteristic (ROC) curves. This curve shows the true-positive rate against the false-positive rate with different classification thresholds. A classification model with AUC = 1 is a perfect classifier, and the classifier with AUC = 0.5 is generally a random selector. RMSD is calculated with respect to the crystal structure.

From Table 1, we observe that MedusaNet performs better than all other approaches in terms of AUC, indicating that MedusaNet has a better pose classification capability. Both CNN models (AtomNet and MedusaNet) can generate poses with less average RMSD compared to the poses proposed by the MedusaDock. The accuracy of MedusaDock is not available since the MedusaDock energy score cannot be used to determine if a pose is good or not. Instead, it only indicates that a pose with a less MedusaDock energy score has potentially less RMSD value than a pose with a higher MedusaDock energy score for the same protein. Although the random forest approach results in good accuracy and average RMSD, it has a roughly 0.5 AUC. In fact, the random forest approach nearly predicted all poses as negative poses and only selected a few as good poses. Since most poses in the data set are bad poses, such selection will get a good accuracy while the AUC is nearly 0.5, which indicates that random forest failed in this task.

In the following sections, we show that not only could our proposed CNN models achieve better scoring capability for protein–ligand complexes, it can further be leveraged to guide the traditional energy-based docking software to be more efficient and effective.

CNN-Accelerated MedusaDock Docking with Improved Accuracy.

In MedusaDock, we perform repeated docking attempts with different random seeds until the minimum energy converges, aiming to find the global minimum energy pose. Specifically, we first perform 64 MedusaDock docking attempts at the outset with different random seeds and then record the pose with the lowest energy. Next, we perform another 64 docking attempts and check if any pose with a lower energy could be found in the 65th to 128th docking attempts. In general, we run twice more attempts of MedusaDock as long as we find a pose with a lower energy, occasionally resulting in thousands of docking attempts. After the minimum energy is converged, we select the k lowest energy poses as the docking results. In this work, we denote k as the number of poses we select for each test. This strategy to find the converged minimum energy is sometimes inefficient because (1) the number of docking attempts will burgeon significantly if the minimum energy dwindles slowly and constantly, and (2) the sampling process could halt at a local minimum, resulting in a low docking accuracy. We utilize the CNN model to circumvent these two difficulties. Specifically, CNN is used to judge if a pose is a good pose or a bad pose. MedusaDock can stop short of new docking attempts if an expected number of good poses have been generated, which will remarkably curtail the total number of docking attempts. On the other hand, CNN can force MedusaDock to explore new ligand poses if the minimum energy converges and the minimum energy pose is judged as a bad pose, thus preventing the docking procedure from being trapped in a local minimum.

We employ two CNN models (AtomNet and our MedusaNet) to guide the MedusaDock and test both models by combining them with MedusaDock through the aforementioned strategy in a data set composed of 100 randomly selected protein–ligand complexes from PDBbind testing set and compare both CNN-augmented versions of MedusaDock (MedusaDock-AtomNet, MedusaDock-MedusaNet) with the original version of MedusaDock. Recall that we randomly choose the training and testing set for 10 times in the cross-validation, the 100 randomly selected complexes are not used for model training in each validation. We observe that the original MedusaDock performs 745 docking attempts averagely among 10 cross-validation sets to achieve a convergent minimum energy (Figure 3). On the other hand, MedusaDock-AtomNet and MedusaDock-MedusaNet only need to perform 370 and 557 docking attempts, respectively, when 128 good poses are selected. If we only select eight good poses, the number of needed docking attempts can slump to 129 and 148, respectively, as compared to the original 745 attempts. Out of the 100 protein–ligand complexes, MedusaDock, MedusaDock-AtomNet, and MedusaDock-MedusaNet can generate near-native ligand pose (RMSD < 2 Å) for 44.6, 41.3, and 49.4 complexes on average, respectively, when eight poses are selected. If we select 128 poses for each complex, MedusaDock, MedusaDock-AtomNet, and MedusaDock-MedusaNet can propose at least one good pose for 67.7, 68.7, and 71.7 complexes on average, respectively. Further, for a large number of selected poses (e.g., 128), CNN-guided MedusaDock might need more attempts compared to the original MedusaDock. For example, on one validation set, MedusaDock-MedusaNet takes 1.27× attempts compared to the original MedusaDock to select 128 good poses for each complex. This result indicates that, for some proteins, performing docking attempts until the energy converges is not a sufficient condition for termination. From Figure 3a, one can also observe that the error bars of the number of docking attempts are very large, which indicates that some proteins will take many attempts of MedusaDock running while other proteins only require a few attempts. This is because the MedusaDock employs a random searching algorithm (Monte Carlo). For some proteins, MedusaDock finds good poses at early attempts, while for other proteins, the good poses are unfortunately proposed in the later attempts. Also, for some complex proteins, MedusaDock can only find a few good poses in total. As the results, we have to run much more attempts for the “Unlucky” proteins than those “lucky” proteins.

Investigating the Effect of Different RMSD Thresholds on Docking.

In the previous evaluation, we use 2 Å as the RMSD threshold to determine if a binding pose is good or not. A pose with an RMSD greater than 2 Å is considered as a bad pose, while the poses with RMSD values less than 2 Å are good poses. However, for different complexes, the RMSD threshold might be different. In this section, we evaluate the performance of our framework with different RMSD thresholds (Figure 4) and check if our approach succeeds with all RMSD thresholds. More specifically, we set the RMSD threshold to 2.00, 2.25, 2.50, 2.75, and 3.00 Å, and perform the convergent docking in the same test set composed of 100 protein–ligand complexes. Again, we applied the cross-validation for 10 times and reported the average, minimum, and maximum results across 10 validation sets for each experimental setting.

In general, for both approaches, the number of successful proteins increases as the RMSD threshold increases. For nearly all RMSD threshold and all combinations of the number of selected poses, MedusaDock-MedusaNet performs better than the original MedusaDock. The number of docking attempts for MedusaDock-MedusaNet increased as the number of selected poses increased. In most cases, MedusaDock-MedusaNet performs less number of attempts compared to the original MedusaDock except on some validation sets. Therefore, our framework achieves significant acceleration while ensuring more protein–ligand complexes find at least one good pose in the designated number of selected poses.