GNINA 1.0: molecular docking with deep learning — TEI-tagged view

Andrew T. McNuttpersName, Paul FrancoeurpersName, Rishal AggarwalpersName, Tomohide MasudapersName, Rocco MelipersName, Matthew RagozapersName, Jocelyn SunseripersName, David Ryan KoespersName (2021date). “GNINA 1.0: molecular docking with deep learning.” Journal of Cheminformaticstitle·j.

Object of study — GNINAterm 1.0 and convolutional-neural-network scoringterm inside a molecular-docking workflow

Research question — whether CNN rescoringterm improves redockingterm and cross-dockingterm pose ranking over Vinaterm while retaining practical runtime and generalizing to unseen protein-ligand systems

Contribution — an open-source GNINAterm 1.0 release, a five-model default CNN ensemble, tuned defaults, and controlled evaluations across four docking tasks

Full text — tagged

GNINAterm 1.0: molecular dockingterm with deep learningterm Abstract Molecular dockingterm computationally predicts the conformationterm of a small molecule when binding to a receptor. Scoring functionsterm are a vital piece of any molecular dockingterm pipeline as they determine the fitness of sampled poses. Here we describe and evaluate the 1.0 release of the Gninaterm docking software, which utilizes an ensemble of convolutional neural networksterm (CNNs) as a scoring functionterm. We also explore an array of parameter values for Gninaterm 1.0 to optimize docking performance and computational cost. Docking performance, as evaluated by the percentage of targets where the top pose is better than 2Å root mean square deviation (Top1term), is compared to AutoDock Vinaterm scoringterm when utilizing explicitly defined binding pockets or whole protein dockingterm. Gninaterm , utilizing a CNN scoringterm function to rescore the output poses, outperforms AutoDock Vinaterm scoringterm on redockingterm and cross-dockingterm tasks when the binding pocket is defined (Top1term increases from 58% to 73% and from 27% to 37%, respectively) and when the whole protein defines the binding pocket (Top1term increases from 31% to 38% and from 12% to 16%, respectively). The derived ensemble of CNNs generalizes to unseen proteins and ligands and produces scores that correlate well with the root mean square deviation to the known binding pose. We provide the 1.0 version of Gninaterm under an open source license for use as a molecular dockingterm tool at https://github.com/gninaterm/gninaterm .

The online version contains supplementary material available at 10.1186/s13321-021-00522-2.

Keywords: Molecular dockingterm, Deep learningterm, Structure-based drug design Introduction Molecular dockingterm is a computational procedure in which the non-covalent bondingterm of molecules, e.g. a protein receptor and a ligand, is predicted. This prediction outputs the conformationterm and, usually, the binding affinityterm of the small molecule in its predicted minimal energy state and is used to virtually screen large libraries of compounds [ 1 – 3 ]. Docking is composed of two main steps: samplingterm and scoringterm. Samplingterm refers to an extensive search of the conformational spaceterm of the molecules being docked. This conformational spaceterm is vast, due in part to both the receptor and ligand being flexible allowing for each molecule to adjust its shape due to the influence of the other. In order to constrain this large conformational spaceterm, the receptor is typically kept rigid. The other vital piece of molecular dockingterm is the scoring functionterm. Every sampled pose is evaluated by the scoring functionterm for its fitness. The fitness determines the conformationsterm that are retained from samplingterm and is used to rank the retained poses in the order of their likelihood of being correct. The final output of docking is a set of ranked poses of the docked molecule.

Determination of the correct binding pose of a small molecule is a prerequisite for determining its binding affinityterm and affords the opportunity to utilize the pose for lead optimization. Correct evaluation of binding affinityterm is critical for downstream tasks such as virtual screening or for determining if a compound is important for more experimental analysis. Molecular dockingterm must compute a pose and a binding affinityterm quickly for it to beneficial when millions of ligands are being queried in a drug discovery pipeline [ 1 ]. Samplingterm the entire conformational spaceterm of a molecule is computationally demanding; therefore, we compromise on the speed and accuracy of docking to provide poses that are close to native while not requiring the full search of conformational spaceterm. This compromise requires docking software to focus on the accuracy and ranking power of the scoring functionterm to highlight low energy conformationsterm and reduce extraneous samplingterm.

Scoring functionsterm provide a mapping from the conformational spaceterm of the ligand and receptor to the set of real numbers so that poses may be ranked. Typically, scoring functionsterm are grouped into three categories; knowledge based, physics based, and empirical [ 1 ]. Knowledge based scoring functionsterm leverage the statistics from a set of structural binding data. A number of geometric properties are computed from structures of protein-ligand complexes, such as atom-atom pairwise contacts. The calculated frequencies can be used in a method such as potential of mean force (PMF), which creates a potential based on the Boltzmann distribution of the properties, to calculate the score of a pose [ 4 , 5 ]. Knowledge based scoring functionsterm can be biased by features present in their training sets though calculations of scores are quick at test time [ 1 ]. They require a large database of known structures and can be difficult to interpret when trying to understand a score [ 6 ]. Physics based scoring functionsterm, often referred to as force fieldsterm, utilize physically derived energetics of interactions to compute scores. The final score is a summation of energy terms such as Coulombic and Van der Waals forces [ 7 ]. Accuracy of physics based scoring functionsterm are limited by their complexity and the assumptions we place on the fundamental forces dictating interactions between molecules, though understanding of these forces is continually increasing [ 8 ].

Empirical scoring functionsterm address the limitations of physics based scoring functionsterm by using a combination of manually selected energy terms. Rather than giving each energy term identical weighting, the weights of each term are determined via a fit to experimental data. A large proportion of docking software use empirical scoring functionsterm, including X-Score, AutoDock Vinaterm, and ChemScore [ 9 – 11 ]. Unlike knowledge based scoring functionsterm, empirical and physics based scoring functionsterm may be easily interpreted to determine the contributing factors of a given score since each energy term can be individually queried. Fitting empirical scoring functionterm·def requires a plethora of experimental structural data and prevents the combination of terms from separately trained scoring functionsterm. The three categories of scoring functionsterm are limited to features extracted from structural information and often assume there is a linear relationship between the features and the binding affinityterm. AutoDock Vinaterm (called “Vinaterm”) utilizes an empirical scoring functionterm·def explicitly tuned to structural data [ 10 ]. The Vinaterm scoring functionterm is a weighted sum of atomic interactions. Steric, hydrophobic, and hydrogen bonding interactions are calculated and are adjusted by the number of rotatable bonds to account for entropic penalties. The weights of the terms were computed via a non-linear fit to structural data. Nguyen et al. [ 12 ] show that Vinaterm can more accurately predict the binding pose than its predecessor, AutoDock 4 [ 13 ]. Vinaterm demonstrates the power of modelling non-linear relationships with its increased docking performance. Therefore, we search for alternative scoring functionsterm that are able to model non-linear relationships between inputs.

Machine learning (ML) represents another growing class of scoring functionsterm [ 8 ]. ML algorithms learn arbitrary relationships between observations and outputs while classical scoring functionsterm assume a specific functional form [ 14 ]. There has been considerable progress in other biomedical fields with the utilization of ML models [ 15 ]. However, machine learning algorithmsterm require a large amount of data to properly generalize to unseen information. The last 20 years has seen a noteworthy increase in the quantity of available protein structures [ 16 ]. A plethora of databases annotate structural data with experimental binding affinityterm data, including PDBbind and BindingDB [ 17 – 19 ]. This information has been utilized to leverage machine learning algorithmsterm as scoring functionsterm. A number of traditional ML approaches have been used as scoring functionsterm, including random forests (RF-Score [ 20 ] and SFCScore [ 21 ]), support vector machines (SVR-Score [ 22 ], ID-Score [ 23 ], SVR-EP [ 24 ]), artificial neural networks (NNscore [ 25 ] and BsN-Score [ 26 ]), and gradient boosted decision trees (BT-dock [ 27 ] and ESPH T-Bind [ 28 ]). These ML methods have been able to match or exceed existing traditional scoring functionsterm. ML methods allow a more robust fit to training data, but are limited to features manually extracted from structural data.

Deep learningterm (DL) methods allow direct inference of features from inputs. They learn a representation of the inputs via layers of simple, non-linear models which transform the representation to higher abstractions to learn complex functions [ 29 ]. DL methods have demonstrated success in a variety of fields, such as computer vision and natural language processing [ 30 , 31 ]. In recent years, there has been significant progress with DL methods in the drug discovery field with many models employing a convolutional framework. Convolutional neural networksterm (CNN) leverage convolutions to infer features directly from input tensors, usually images. CNNs have shown potential in virtual screening (AtomNet [ 32 ], DeepVS [ 33 ], RagozapersName et al. [ 34 ]) and binding affinityterm prediction (PotentialNet [ 35 ], K DEEP [ 36 ], Pafnucy [ 37 ]). A number of methods have been proposed to capitalize on the power of DL scoring functionsterm. MedusaNet uses a CNN within the docking pipeline to guide the samplingterm of the base docking method [ 38 ]. The base docking method, Medusa, provides a variety of ligand poses. The CNN evaluates the 3D coordinate representation of the poses to determine if a pose should be retained. Nguyen et al. [ 39 ] describe a generative adversarial network (GAN) for pose prediction. Their network utilizes an encoder with low-dimensional mathematical representations of the protein-ligand complex and a decoder utilizing convolutional layers to generate and rank ligand poses for the D3R grand challenge. MasudapersName et al. [ 40 ] use a receptor structure as the prior to their GAN to sample novel ligands appropriate to the identified binding siteterm.

Previous work has largely evaluated deep learningterm protein-ligand scoringterm on already generated poses. In this work, we describe and comprehensively evaluate version 1.0 of the Gninaterm molecular dockingterm software, a fork of Sminaterm [ 41 ] and AutoDock Vinaterm [ 10 ] that supports CNN scoringterm as an integral part of the docking workflow. Gninaterm is evaluated here for its ability to properly score and rank binding poses for protein-ligand complexes. We describe how the default settings which balance docking accuracy and runtime were determined, including the selection of a default ensembleterm of CNN models. Performance of Gninaterm is evaluated for the redockingterm, cross-dockingterm, flexible dockingterm, and whole protein dockingterm tasks and is found to significantly outperform Sminaterm /Vinaterm in all cases.

Methods The docking pipeline of Gninaterm (Fig. 1 ) is described in detail, providing background for the derivation of default usage. A default CNN ensemble is selected for optimizing the docking performance and runtime of the docking pipeline. This ensemble is then used to investigate the different CNN scoringterm options available to the user, followed by a thorough investigation of the docking parameters. Finally, we examine both the generalizability and scoringterm power of Gninaterm .

The Gninaterm samplingterm and scoringterm algorithm shown with relevant commandline parameters and the scope of CNN scoringterm Molecular dockingterm pipeline Gninaterm is a fork of Sminaterm [ 41 ] which is a fork of Vinaterm[ 10 ]. The docking pipeline of Gninaterm utilizes the enhanced support for scoringterm enabled in Sminaterm to allow the use of CNNs as scoring functionsterm. In typical usage, Gninaterm is provided with a receptor structure, a ligand structure, and a specification for a binding siteterm on the receptor (Fig. 2 ).

Example Gninaterm usage Open Babel [ 42 , 43 ], a chemical toolbox allowing the reading and writing of over 100 chemical file formats, is used for parsing the inputs, allowing commonly used structural data formats (e.g. PDB, sdf, mol, etc.) as well as gzipped versions of such files to be used. The binding siteterm can be specified as a Cartesian box or by providing a ligand file (autobox_ligandterm). When autobox_ligandterm is used to define the binding siteterm, a rectangular prism is constructed using the minimum and maximum values for the x , y , and z coordinates of the ligand to which additional spacing (autobox_addterm) is added in every dimension (Figs. 1 , 2 ). In Gninaterm , if any side of this auto-generated box is smaller than the longest distance between any two atoms in the ligand, then those sides are extended to that longest distance, ensuring that the ligand can rotate freely within the defined box without incurring an out-of-box penalty that is applied to all docked poses to constrain them to the specified binding siteterm search space.

Next, the scoring functionsterm are setup. Similar to Sminaterm , if the user opts to not use the CNN scoringterm function for a part of the pipeline they can specify their own empirical scoring functionterm·def to Gninaterm or choose one of the built-in scoring functionsterm, i.e. Vinaterm, Vinardo [ 44 ], etc. CNN scoringterm functions can be specified by providing model and/or weights files or by selecting a built-in model. The available built-in CNN models include crossdock_default2018, dense, general_default2018, redock_default2018, and default2017, each of which is trained using different training data and/or a different model architecture. Additionally, each model, except for default2017, has five variants that are trained on the same data and have the same architecture, but are initialized with a different random seed. For each CNN model type, we refer to the variant with the highest docking performance as the base model name, and the remaining variants are given sequential numbers (i.e. general_default2018, general_default2018_1, general_default2018_2, general_default2018_3, and general_default2018_4); the ensemble of these five variants is denoted with ‘_ensemble’ (i.e. general_default2018_ensemble). The architecture and training of these models are described elsewhere [ 34 , 45 ]. CNN calculations are performed using the cuDNN accelerated Caffe deep learningterm framework [ 46 ]. These models are trained to predict both a pose scoreterm (a probability that the pose has a low root mean square deviation (RMSDterm) to the binding pose) and the binding affinityterm (pK). The pose scoreterm is used for all pose optimization tasks. The scoring functionterm used in each step of the Gninaterm pipeline is determined by cnn_scoringterm (defined below), defaulting to only using the CNN scoringterm for the final resorting of ligand conformationsterm and the empirical scoring functionterm·def everywhere else in the pipeline.

The docking procedure uses Monte Carlo samplingterm to search the ligand conformational spaceterm. exhaustivenessterm (default 8) defines the number of Monte Carlo chains that are run for the ligand. Chains are run in parallel using up to cpu threads. The number of steps for the Monte Carlo chains are calculated based on the number of mobile atoms and the number of degrees of freedom within the ligand. This calculation can be ignored and the user may explicitly specify the number of steps for the Monte Carlo chains with num_mc_steps. Each step of a Monte Carlo chain mutates the ligand by randomly selecting one of the following operations: random translation, random rotation of the entire molecule, or randomly setting the torsional angles of the ligand. The Monte Carlo process selects random resetting of the torsional angles with a higher probability than the other mutations. After the mutation, an approximate energy minimization of the ligand is performed. If an empirical scoring functionterm·def is specified for guiding samplingterm, this minimization is performed using a fast, grid-based approximation of that scoring functionterm. A grid is pre-calculated for each ligand atom type and a single atom of that type is used to calculate values for each grid point. The full ligand is scored by interpolating values from the grid for each of its atoms and summing the result. If a CNN scoringterm function is specified, no such approximation is used since CNN scoringterm is not additive with respect to the individual atoms [ 47 ]. The score for the minimized conformationterm determines if it will be accepted using the Metropolis acceptance criterionterm. Each Monte Carlo chain retains its top scoringterm ligand conformationsterm, and the number retained is user configurable with num_mc_savedterm (default 50).

Following the completion of Monte Carlo samplingterm, the saved conformationsterm from each Monte Carlo chain are aggregated and the num_mc_savedterm top scoringterm conformationsterm are retained for further analysis. The top scoringterm conformationsterm are then refined. If an empirical scoring functionterm·def is specified for refinementterm (cnn_scoringterm set to “none” or “rescore”), a functional formulation, rather than a grid approximation, of the scoring functionterm is used to carry out refinementterm in order to get the locally optimal pose. Refinementterm shifts the ligand pose to a local energy minimum using the gradients of the scoring functionterm. After the ligand pose has been refined, the final affinities and scores are calculated for the pose using the specified CNN models and/or the specified scoring functionterm. Finally, the top scoringterm conformationsterm are sorted by pose_sort_order (default “CNNscoreterm”) and output with Open Babel in the user-specified format if an output file was provided.

The usage of the CNN models within the docking pipeline can be selected by the user using cnn_scoringterm (Fig. 1 , 2 ). If “none” is selected for cnn_scoringterm, then the CNN models are not used at all in the docking pipeline, making the pipeline essentially identical to the Sminaterm pipeline. The only differences are that Sminaterm computes with double (64 bit) precision rather than the single (32 bit) precision used by Gninaterm and Sminaterm does not include the autobox_extend parameter when creating the samplingterm box. When “rescoringterm” is selected (the default), the specified CNN models are not used until the final sorting of resultant ligand conformationsterm after their refinementterm with the empirical scoring functionterm·def. In this case, the specified CNN models are used to score each of the ligand conformationsterm and the output ligand conformationsterm are resorted based on the score calculated by the CNN model(s). The “refinementterm” option utilizes the CNN for the refinementterm of the ligand poses after they have been selected by the Monte Carlo chains and then sorts the refined ligand conformationsterm by the CNN score for output. Finally, the “all” option utilizes the CNN for all aspects of the docking pipeline including the minimization within the Monte Carlo chains, the refinementterm after the Monte Carlo chains, and the sorting of the final output.

Flexible dockingterm Molecular dockingterm is often performed using a rigid protein targetterm and only samples the conformational spaceterm of the ligand, as described above. This is a good approximation when redockingterm to a receptor that does not undergo conformational change upon ligand binding, but it fails to accurately represent the biochemical aspects of the system when the protein undergoes significant structural changes upon binding [ 48 ]. Allowing the whole target to be flexible is too computationally expensive for docking; however, some protein flexibility can be accounted for by samplingterm the conformational spaceterm of side chainsterm in the binding siteterm [ 49 ] while leaving the backbone fixed. The local flexibility of the binding siteterm could be especially beneficial for cross-dockingterm, i.e. docking a ligand to a non-cognate receptor. Assuming a rigid receptor during cross-dockingterm is not realistic when potentially very different ligands are docked to the same receptor.

Gninaterm allows the samplingterm of side chainterm, but not backbone, conformational spaceterm. More specifically, side chainterm torsion angles are sampled as part of the Monte Carlo search. Side chainterm flexibility can be specified manually or semi-automatically in several different ways: flexible side chainsterm can be defined in a PDBQT file (flex parameter), they can be selected using a comma-separated list of residue identifiers (chain, residue number, and, optionally, insertion code; flexres parameter), or they can be selected based on their distance from a given ligand (flexdist and flexdist_ligand parameters).

The last option, used in this work, is similar to the autobox_ligandterm option described above. The flexdist parameter allows the specification of a threshold distance from the flexdist_ligand. If a residue has any side chainterm atoms that are within this distance of the specified ligand, then the entire residue side chainterm is marked as flexible.

Given the increased computational cost of samplingterm side chainterm conformationsterm, two additional options are provided for flexible dockingterm. flex_limit allows users to specify a hard upper bound for the number of flexible side chainsterm. If flex_limit is exceeded for a particular system, docking does not run and Gninaterm terminates with a warning. This option is particularly useful to avoid possible bottlenecks in large virtual screening tasks. The less strict flex_max option allows users to prune the list of flexible side chainsterm in order to make the calculation more manageable for large systems; if there are more than flex_max side chainsterm identified as flexible, only the conformational spaceterm of the flex_max closest to flexdist_ligand is sampled during docking.

Flexible side chainsterm are selected at the very beginning of the docking procedure and the selection is not updated during samplingterm. When using autobox_ligandterm to automatically define the search space as described above, flexible side chainsterm are included in the calculation of the box bounds. The usage of CNN models for docking with flexible side chainsterm can be selected by the user in the same way it is done for docking with a rigid receptor, as described above.

Data There are two primary ways to evaluate molecular dockingterm: redockingterm a cognate ligand to its receptor and docking a ligand to a non-cognate receptor (cross-dockingterm). In order to best evaluate the performance of Gninaterm for molecular dockingterm, we evaluate its performance on both of these tasks. Redockingterm the cognate ligand demonstrates the samplingterm and scoringterm power of the molecular dockingterm pipeline, as the RMSDterm from the crystal pose can be measured to exactly determine the accuracy of the produced poses. Analysis of redockingterm requires a set of high quality structures in which the native binding pose of the ligand has been solved. For this purpose we utilize the PDBbind refined set v.2019term [ 18 ]. The PDBbind database is a curated set of protein-ligand complexes containing both structural information and binding affinityterm. The PDBbind database is updated annually with new experimentally determined structures annotated with binding affinityterm data. The refined set is a subset of the entire PDBbind dataset that retains only the structures with resolution higher than 2.5 Å, high quality affinity measurements, and binary protein-ligand complexes. The 2019 release of the refined set contains 4,852 high quality crystal structures of native protein-ligand binding poses.

However, redockingterm is not the normal use case of a molecular dockingterm pipeline. Docking will ordinarily be performed on proteins-ligand pairs that have no co-crystallized structure. Often a new ligand will be docked into a receptor whose co-crystallized ligand is a different molecule. Wierbowski et al. [ 50 ] recently published a dataset that provides a benchmark precisely for this task. This cross-docking datasetterm provides a meaningful method for the evaluation of the ligand RMSDterm from the known and predicted poses. A reference structure is selected for the protein, then the “known” binding pose of a ligand is defined by the ligand’s position when an alignment is performed between the reference structure and the ligand’s co-crystal receptor structure. The dataset is composed of 94 unique protein binding pockets and 4,399 unique ligands, with an average of 46 ligands per target.

Both of the datasets were filtered to ensure the protein-ligand structures can be parsed by Gninaterm . ProDy [ 51 ] was used to separate the complexes into protein and ligand files while removing any water or other extra crystallized molecules. Our goal is the binding pose prediction of a small molecule at its target site, therefore we utilized RDKit [ 52 ] to filter both the redockingterm and cross-docking datasetsterm to include only ligands with molecular masses greater than 150 Da and less than 1000 Da. Any ligand that was not able to be parsed with RDKit was also removed. In addition, we visually inspected the ligands of the cross-docking datasetterm with Pymol [ 53 ] to ensure that each ligand for a given binding pocket overlapped. If a ligand was non-overlapping it was marked as problematic. Lastly, for the UROK pocket, 4ZKO was removed as the receptor did not align properly and the ligands from 4MNW and 4MNX were also marked as problematic due to their binding pose depending on a chain from the kinase being in a protein-protein complex rather than a holo structure as in the other members of the pocket. The final filtered datasets were composed of 4,260 and 820,280 protein-ligand pairs for the redockingterm (PDBbind [ 18 ]) and cross-dockingterm (Wierbowski et al. [ 50 ]) datasets, respectively.

Due to the combinatorial size of the cross-docking datasetterm, more filtering was required to make computational time tractable. Each unique protein binding pocket can be used to group the protein-ligand complexes as all of the receptor structures in the group share a common pocket. These groupings allow the cross-docking datasetterm to be downsampled for each pocket. Each pocket was either kept in full, or reduced to a random sample of 100 receptor-ligand pairs, whichever was smaller. We then removed any of the ligands that were previously listed as problematic. The downsampling results in 7970 protein-ligand pairs and 92 unique protein binding pockets in the cross-docking datasetterm, where in no case is a ligand paired to its cognate receptor. Notably, 2 pockets were removed (ALDR and CP3A4) due to all of the ligands either failing the filtering, or being marked as problematic. We evaluate the docking performance differences for various fractions of the total complexes per pocket for both Gninaterm and Sminaterm to ensure the downsampling does not bias the performance of either software (Additional file 1 : Fig. S1).

Evaluation metrics TopNterm : We evaluate docking performance by examining the output poses. The Open Babel obrms tool [ 42 , 43 ] is used to determine the RMSDterm from each output pose to the binding pose. In redockingterm, the binding pose is defined by the crystal pose. In cross-dockingterm, we define the binding pose as the conformationterm of the ligand when the cognate protein structure is aligned with the reference structure. If the RMSDterm to the known binding pose is less than 2 Å then we consider the pose to be “good.” The percentage of systems with a good pose ranked in the top N (TopNterm) is reported for redockingterm. In the case of cross-dockingterm, we consider all poses of a ligand across an ensemble of non-cognate receptors. We calculate the TopNterm for each target individually and report the average across the 92 systems. Averaging is performed to avoid over-representing targets with a large number of ligands. This metric can be computed for any number of output poses, computing only the top pose is Top1term.

Avg Time Per Systemterm·var: Properly benchmarking the time for the various CNN scoringterm options on both the redockingterm and cross-docking datasetsterm would take a significant amount of time, so we utilize a filtered version of the PDBbind v.2016 [ 54 ] core set to determine runtime. This set consists of 263 complexes. The average runtime is calculated per system using the hyperfine benchmarking tool [ 55 ], computing the runtime for a minimum of 5 docking runs for each system (additional runs are automatically performed by hyperfine if high variability is encountered). We then average docking runtime over the whole PDBbind core set, to provide the average time for one docking run of the PDBbind core set. Timing evaluations were done on a dual 16-core 2.3 GHz Intel Xeon 5218 with 96 GB of RAM and a 11 GB RTX 2080 Ti GPU using CUDA 10.2 and cuDNN 7.6.5. All benchmarking runs were done on an otherwise unloaded system with four cores requested (cpu=4).

Default model selection When using CNN scoringterm, the user can utilize a single CNN model or an ensemble of CNN models where the final score is an average of each CNN model’s score. It has been shown that ensembles of predictors improve performance over a single predictive model [ 56 ]. However, due to the high computational cost of applying CNN models in comparison to empirical scoring functionsterm, it is desirable to select a subset of the available models that provides improved scoringterm while limiting computational cost. The default CNN model ensemble was selected using a greedy forward algorithm. The ensemble was built in an iterative process using the “rescoringterm” option for CNN scoringterm so as to minimize the computational time for each iteration. In each round of the selection, models were chosen for their Top1term performance. We evaluate all of the versions of each CNN model (i.e. dense, dense_1, dense_2, dense_3, dense_4). In the first round of selection, all of the CNN models were tested for their individual ability to predict Top1term on both the cross-dockingterm and redockingterm datasets. The next round of selection required testing of all two-model combinations with the model selected in the first step. Model selection continued, exploring all possible combinations of the built-in CNN models. The selection process concluded after five CNN models were selected for inclusion in the default ensembleterm.

Default CNN scoringterm method Gninaterm allows the usage of CNN scoringterm in various steps of the molecular dockingterm process (Fig. 1 ). The CNN scoringterm option allows the user to change how the CNN is used to evaluate a ligand pose. If the CNN is not used at all in the scoringterm process (“none” option), then the molecular dockingterm pipeline is essentially the same as Sminaterm . The “rescoringterm” option has the lowest computational cost of the options that utilize the CNN models. With this option CNN models are used to score and re-sort the ligand conformationsterm selected and refined by the non-CNN scoring functionterm, defaulting to the Vinaterm scoring functionterm. For additional computational cost, the “refinementterm” option can be specified to use the CNN models for the refinementterm of the ligands after they have been selected by the Monte Carlo chains. In addition to refining the ligand conformationsterm, the CNN models are used to score and resort the output poses, as the “rescore” option does. However, the Monte Carlo chains continue utilizing the non-CNN scoring functionterm. The “all” option uses the CNN models as the scoring functionterm throughout the course of the molecular dockingterm procedure and has the highest computational cost by orders of magnitude. The CNN model is used for the selection process within the Monte Carlo chains, the refinementterm process after Monte Carlo selection, and the scoringterm and resorting of the poses before output. This option is very computationally intensive as the CNN is regularly queried for the energy of a particular conformationterm during the Monte Carlo samplingterm procedure.

In addition to the values allowed for cnn_scoringterm, the user is provided with cnn_empirical_weight to combine the non-CNN and the CNN scoringterm functions. Using mix_emp_force the refinementterm of the ligand poses can be computed with a linear combination of the CNN gradients and the non-CNN scoring functionsterm force. mix_emp_energy uses the same linear combination of the scoring functionsterm for computing the score of a given pose. The weighting of the Vinaterm scoring functionterm within the linear combination is selected with cnn_empirical_weight (default 1.0).

The default usage of the CNNs within the Gninaterm docking pipeline demands high accuracy while limiting computational costs. Therefore, each of the CNN scoringterm options is evaluated for both docking performance and runtime. Docking performance can be calculated on both the redockingterm and cross-docking datasetsterm via TopNterm. To this end, Gninaterm is used with one of the CNN scoringterm options with the default CNN ensemble selected above to compute 9 ligand poses for each protein ligand system. We also investigate various values of cnn_empirical_weight while using the “refinementterm” option to determine if a linear combination of the non-CNN and CNN scoringterm functions provides greater docking performance. We use both mix_emp_energy and mix_emp_force with the cnn_empirical_weight parameter to use the linear combination of Vinaterm and the default CNN ensemble for both refinementterm of poses and final scoringterm.

Parameter exploration Gninaterm has many parameters that alter the molecular dockingterm pipeline (Fig. 2 ). A default value is found for each of these parameters to provide the best all around default behavior. In exploring the various values of the parameters we find the optimal value for both the redockingterm and the cross-docking datasetsterm. Additionally, we also consider that the user may not know the exact location of the binding pocket on the receptor and may have to use the entire protein as the autobox for the docking pipeline. Therefore, our setting exploration considers the case in which the specific binding pocket is known and the case in which the whole protein is used to define a binding box. Some parameters directly impact the samplingterm procedure, such as exhaustivenessterm, autobox_addterm, num_modesterm, num_mc_savedterm, and min_rmsd_filterterm·var which are described below. These parameters control, to some extent, the extensiveness of the search during the Monte Carlo samplingterm procedure. As previously stated, exhaustivenessterm determines the number of Monte Carlo chains run during the samplingterm procedure. autobox_addterm increases the size of the binding box that the Monte Carlo chains sample. num_modesterm determines the number of ligand poses output by Gninaterm at the completion of the docking procedure. This is separate from num_mc_savedterm which defines the number of ligand poses saved for each Monte Carlo chain. The number of ligand poses retained after all of the Monte Carlo chains are completed is determined by either the number of modes or the number of Monte Carlo saved, whichever is larger. After all of the Monte Carlo chains have completed and the poses have been refined and sorted, the RMSDterm between all pairs of ligand poses is calculated. min_rmsd_filterterm·var removes one pose from a pair if the RMSDterm of the pair is less than the value of the parameter. This ensures the poses returned by the docking procedure are all different from one another. When using the CNN scoringterm function, another setting is how many different rotations of the protein-ligand complex the CNN is able to see for each conformationterm, cnn_rotationterm·var.

Evaluations for all of the parameters are carried out on both the redockingterm and cross-docking datasetsterm. Each setting is varied individually using the default CNN ensemble determined above. The values explored for each parameter are defined in Tables 1 and 2 . Values were explored around the previously set defaults for each of the parameters. Each value produced a set of poses, which was used to calculate a TopNterm.

Parameters explored when the binding pocket has been defined Parameters explored when the binding pocket is not known and the whole protein is used for docking CNN scoringterm performance All of the CNN models were trained on some subset of the cross-dockingterm and redockingterm datasets. Generalization can be evaluated by determining the performance of the CNN scoringterm functions on the subset of the datasets that were not seen in training. This evaluation is carried out by removing any protein or ligand contained within the training data of the CNN models. All of the CNN models have been trained on different sets, so fully testing generality requires the removal of all of the proteins and ligands in the training set from both the redockingterm and cross-docking datasetsterm. We removed all PDB IDs that were in the training sets of the CNN models [ 18 , 45 ], leaving 441 and 178 protein-ligand pairs for the redockingterm and cross-dockingterm sets, respectively.

The CNN models output both a CNNscoreterm and a CNNaffinityterm for each of the conformationsterm output by Gninaterm . CNNaffinityterm is the affinity of the docked complex as determined by the CNN, this metric has been evaluated in a previous work [ 45 ]. The CNNscoreterm is a value between 0 and 1 that is used to rank the poses of the ligand, where a score of 1 denotes a perfect ligand pose. We would like to investigate if there is a correlation between high scores and low RMSDterm to the crystal pose. For ease of analysis, we only consider the top ranked pose. Using the top ranked pose for each complex, we investigate how filtering the poses by their CNNscoreterm can affect the percentage of poses in which the RMSDterm to the crystal pose is less than 2 Å.

Results Sminaterm Comparison Gninaterm is a fork of Sminaterm that allows the utilization of CNN models as scoring functionsterm. Therefore, without the use of the CNN models Gninaterm should function exactly as Sminaterm does. However, unlike Sminaterm , Gninaterm does computation with single (32 bit) precision rather than double (64 bit) precision due to the need to shift calculations to the GPU for efficient CNN scoringterm. Therefore, we ensure that the use of single precision does not negatively affect the docking power of the pipeline. The effect of this precision change can be evaluated by running Gninaterm without using any CNN scoringterm and with autobox_extend turned off, allowing us to compare to Sminaterm docking results. We consider only redockingterm results as cross-dockingterm results would require significantly more computation, and identification of differences due to precision can be done using only redockingterm results. Results for redockingterm do not show a significant difference for the output poses. A majority of the output poses are exactly the same, with slight differences seen for some output poses (Additional file 1 : Figs. S2 and S3).

Default model selection The iterative process for construction of the default CNN ensemble, denoted Default Ensembleterm, is shown in Additional file 1 : Table S1. The five selected models are dense, general_default2018_3, dense_3, crossdock_default2018, and redock_default2018. We now evaluate the docking performance boost that this ensemble provides over any single CNN model type (e.g. crossdock_default2018, dense, etc.) or an ensemble of the same CNN model type (e.g. crossdock_default2018_ensemble, dense_ensemble, etc.). We compare the CNN model(s) docking performance by evaluating TopNterm on both the redockingterm and the cross-dockingterm tasks. The CNN models are used in the “rescoringterm” option for the CNN scoringterm to output 9 ligand conformationsterm.

The docking performance of the Default Ensembleterm is compared to the single model options in Fig. 3 . While nearly all models are able to outperform Vinaterm, we can see that the newly selected Default Ensembleterm significantly outperforms all of the single models on both the redockingterm task and the cross-dockingterm task.

Docking using the single CNN models and the newly selected Default Ensembleterm for rescoringterm the output poses. The binding pocket is defined by the known binding ligand. TopNterm is the percentage of targets ranked above or at N with a RMSDterm less than 2 Å When the Default Ensembleterm is compared to the ensembles of each of the individual CNN model types, we see that the Default Ensembleterm is able to outperform all of the ensembles composed of one model type (Fig. 4 ). The ensemble selection procedure determined five CNN models whose combined performance on ranking low RMSD posesterm first beats the performance of the ensemble utilizing all of the built-in models while significantly reducing computational cost. The All Ensembleterm is the ensemble composed of all 16 CNN models built-in to the Gninaterm software. The Default Ensembleterm is able to meet the docking performance of this large ensemble, even when considering cross-dockingterm results per pocket (Additional file 1 : Fig. S4), while only being composed of five models. Reducing the number of models in the ensemble enables the computations to be several seconds faster for an average docking run (Fig. 5 ). This reduction is likely due to the inclusion of only two of the dense models which take the longest to run because of the high number of parameters in the models. The computational speedup afforded by the Default Ensembleterm over the All Ensembleterm increases when no GPU is used for docking (Additional file 1 : Table S2). The computational speed boost can have a significant impact when performing a large number of docking runs or when there is no GPU available for enhanced parallelism of the scoringterm computation.

Docking using the ensemble of each type of CNN model, the full ensemble of CNN models, and the newly selected Default Ensembleterm for rescoringterm the output poses. The binding pocket is defined by the known binding ligand. TopNterm is the percentage of targets ranked above or at N with a RMSDterm less than 2 Å Evaluation of the average time to dock one protein-ligand system from the PDBbind core set v.2016. Top1term is the percentage of top ranked targets with a RMSDterm less than 2 Å When comparing Figs. 3 and 4 we can see that the ensembles composed of the individual model types are able to outperform their single model counterparts. We therefore omit the single models for the remaining evaluations.

Default CNN scoringterm method We evaluate the performance of the Default Ensembleterm with the “rescore,” and “refinementterm” options of cnn_scoringterm. The usage of the “all” option was unable to complete on the PDBbind Core set in a reasonable amount of time, so it was not considered. The “all” option requires usage of the CNN scoringterm function for every mutation of the ligand within every Monte Carlo chain with each usage of the CNN scoringterm function having a high computational cost. This leaves us with the “refinementterm” and “rescoringterm” options of cnn_scoringterm. Figure 6 shows the Default Ensembleterm performs nearly as well with either option. We can also see that using cnn_emp_weight with both mix_emp_energy and mix_emp_force does not significantly alter the docking performance when using the “refinementterm” option (Additional file 1 : Figure S5).

Comparing the Default CNN Ensemble for use in only rescoringterm of the poses output by the Monte Carlo chains or the refinementterm of the poses followed by a rescoringterm of the poses. The “refine” option has nearly the same docking performance as the “rescore” option when cross-dockingterm. TopNterm is the percentage of targets ranked above or at N with a RMSDterm less than 2 Å However, from looking at the average time to perform molecular dockingterm for one system we see that “refinementterm” takes an order of magnitude longer than “rescoringterm” (Additional file 1 : Figure S6). Time for performing “rescoringterm” on an average system is similar to the time to perform docking with the Vinaterm scoring functionterm. We find it reasonable to use “rescore” as the default option for the CNN scoringterm due to its docking performance and runtime.

Parameter exploration Changes in the exhaustivenessterm alter the amount of samplingterm that occurs. When the exhaustivenessterm is increased, we see an increase in the performance of docking (Fig. 7 ). This is as expected as more Monte Carlo chains randomly mutating the ligand conformationterm provides the docking procedure with more opportunities to randomly sample the correct pose. However, there are no significant performance gains after a value of 8. An exhaustivenessterm of 16 provides some performance boost, but this boost may be accompanied by a doubling of the computational time. The Monte Carlo chains are evaluated in parallel, but parallelism is limited by the number of cores available to Gninaterm . If the exhaustivenessterm is greater than the CPUs provided to Gninaterm , the number of simultaneously running Monte Carlo chains is equal to the number of CPUs. Therefore, in typical usage of Gninaterm , an exhaustivenessterm level of 8 is sufficient for adequate performance levels while minimizing the computational load when targeting a specific binding siteterm.

Evaluating the role of exhaustivenessterm in the performance of docking with the Default CNN Ensemble by analyzing TopNterm. TopNterm is the percentage of targets ranked above or at N with a RMSDterm less than 2 Å Increasing the num_mc_savedterm parameter increases the chances of samplingterm accurate docking poses, as increasing the number of output conformationsterm from each Monte Carlo chain increases the likelihood of finding the correct pose. However, this will also increase computational time due to the fact that more poses require refinementterm and final scoringterm. As num_mc_savedterm gets closer to 100, we see that the docking performance boost is reduced (Fig. 8 ). Therefore, we select a new default value of 50 for num_mc_savedterm to minimize the computational overhead while still increasing the performance of the docking routine. The num_mc_savedterm and num_modesterm affect one another; the number of poses saved from each Monte Carlo chain is the maximum of the two values. When looking at the first 9 poses, we see an increase in the docking performance with a substantially greater value for the number of output poses (num_modesterm) (Fig. 9 ). This is due to num_modesterm forcing each Monte Carlo chain to output a number of poses greater than num_mc_savedterm. Increasing the default value of num_modesterm to a value higher than 50 (the default value for num_mc_savedterm) will again increase computational overhead, so the default value is set to 9.

Evaluation of the Number of Monte Carlo Saved in the performance of docking with the Default CNN Ensemble by analyzing TopNterm. TopNterm is the percentage of targets ranked above or at N with a RMSDterm less than 2 Å Evaluating a much greater number of modes on the performance of docking with the Default CNN Ensemble by analyzing TopNterm. TopNterm is the percentage of targets ranked above or at N with a RMSDterm less than 2 Å We next evaluate the size of the targeted binding siteterm (box). The autobox_addterm parameter increases the search space for the docking program to be larger than the rectangular prism defined by the autobox_ligandterm input to Gninaterm . In redockingterm (Additional file 1 : Figure S7a), the expansion of the search space decreases performance as the expanded search space increases the potential conformational spaceterm of the ligand. A higher value of autobox_addterm is necessary to ensure that the correct binding siteterm, which may differ from the binding siteterm of the autobox_ligandterm, is included in the search space. This is shown during cross-dockingterm (Additional file 1 : Figure S7b), where only the binding siteterm of the cognate ligand is known. A low value of autobox_addterm artificially improves docking performance when the binding pose is known by unrealistically constraining the search space. We therefore select a default value of 4 for autobox_addterm to keep the box small while still providing room for error in the selection of the binding siteterm.

Changes to the CNN rotations do not significantly change the scoringterm of the Default Ensembleterm (Additional file 1 : Figure S8). The CNN Ensemble is able to determine the correct score for the ligand pose regardless of the rotation of the ligand and protein complex. Altering the value of the minimum RMSDterm filter does not change the results of the docking (Additional file 1 : Figure S9). Filtering out poses with similar conformationsterm increases the diversity of poses that the CNN ensemble ranks. However, the CNN ensemble is able to accurately rank the poses it sees, providing high scores to poses with low RMSDterm to the known binding pose.

Whole protein dockingterm Next, we evaluate the performance of docking when using the whole protein as the defined binding siteterm. Whole protein dockingterm can be used for new protein targets when the true binding siteterm is unknown. The performance of docking is expected to be reduced as the samplingterm space has been significantly increased. When using the whole protein for the samplingterm space, the ligand tends to get stuck at local energy minima which are distant from the actual binding siteterm. Most of the protein surface is not hospitable for ligand binding, so once a potential pocket is discovered, exploring more of the protein surface has a low probability during the Monte Carlo samplingterm procedure. Comparing Figs. 10 to 7 shows that the docking performances of both Vinaterm and the Default CNN Ensemble are reduced from when the binding pocket was explicitly defined. Redockingterm performance decreases from 73% to 38% at Top1term with exhaustivenessterm at 8 while cross-dockingterm decreases from 37% to 16% at Top1term with exhaustivenessterm at 8 (cross-dockingterm per pocket Top1term is shown in Additional file 1 : Figure S10). The larger potential docking space requires more samplingterm to find a ligand conformationterm with low RMSDterm to the known binding pose. As expected, with whole protein dockingterm we see greater increases in performance with increased samplingterm (exhaustivenessterm). By again comparing Figs. 10 with 7 we see that this boost in performance is larger than when the binding pocket is defined explicitly and does not exhibit the same diminishing returns. When increasing exhaustivenessterm from 8 to 16 in whole protein dockingterm Top9 increases from 48% to 58% and 22% to 29% for redockingterm and cross-dockingterm, respectively. If the binding pocket is explicitly defined, the same change in exhaustivenessterm only increases Top9 from 87% to 88% and 53% to 55% for redockingterm and cross-dockingterm, respectively. For this reason, when performing whole protein dockingterm, we recommend setting the value of exhaustivenessterm as high as possible given time constraints.

Increasing the exhaustivenessterm when using the whole protein as the binding box. TopNterm is the percentage of targets ranked above or at N with a RMSDterm less than 2 Å Flexible dockingterm We now test the performance of Gninaterm for docking with flexible side chainsterm. Since receptor flexibility is only useful in the context of cross-dockingterm, we limit our tests to the cross-docking datasetterm; allowing side-chain re-arrangements would only deteriorate the performance in redockingterm given that the receptor is already in the correct conformationterm.

Flexible dockingterm is computationally more expensive than docking with a rigid receptor because of the larger number of degrees of freedom to be sampled. For this reason, we only use the default parameters carefully selected above. In order to define the side chainsterm to be treated as flexible, autobox_ligandterm is also used as flexdist_ligand and flexdist is set to 3.5 Å, which gives a reasonable representation of the protein-ligand binding siteterm [ 57 ]. Therefore, conformationsterm for all side chainsterm with at least one atom within 3.5 Å from flexdist_ligand are sampled during docking.

Out of the 7970 systems in the cross-docking datasetterm, flexible side chainsterm are not identified for 24 protein-ligand complexes (see Additional file 1 : Table T3 for details). Such systems are discarded from the following analysis of flexible dockingterm since they are equivalent to rigid docking. Additionally, 9 systems were problematic when computing RMSDsterm for the flexible residues because of connectivity issues. Disulfide bonds between cysteine residues are allowed to break during samplingterm (with a software warning), which results in a different connectivity in the output file for four systems. In the other five systems, the connectivity between flexible residues in the input and output file was found to be different and therefore they were discarded (see Additional file 1 : Table T3 for details). All nine systems with connectivity issues were removed from the analysis, resulting in a dataset of 7937 cross-docked complexes.

In order to assess the performance of flexible dockingterm compared to rigid docking, we look at RMSDterm differences between ligands docked with both methods in relation to the similarity between binding pockets of the cognate and docking target receptors. The similarity of the binding pockets is assessed via side chainsterm RMSDterm between the docking target (non-cognate receptor) and the cognate receptor, which we denote target-cognate RMSDterm. This is distinct form either the side chainterm RMSDterm between docking input and output (target-pose RMSDterm) and the side chainterm RMSDterm between cognate receptor and docking output (cognate-pose RMSDterm). The target-cognate RMSDterm is computed by finding the best match between residues in the target and cognate receptors that are within 3.5 Å from the ligand being docked.

Figure 11 a shows the difference in RMSDsterm for the ligand top pose between flexible and rigid docking versus the target-cognate RMSDterm. For higher target-cognate RMSDsterm, indicating differences in the binding pocket between the target and cognate receptor, one would expect flexible dockingterm to perform better. However, as we can see, the difference in ligand RMSDterm between flexible and rigid docking for the top pose varies widely between systems and there in no clear advantage in flexible dockingterm. The overall RMSDterm distributions for the ligand top poses (Additional file 1 : Figure S11) are fairly similar, with slightly more systems with low RMSDterm for rigid docking than flexible dockingterm.

Comparison between rigid and flexible dockingterm with the default Gninaterm parameters: ( a ) ligand RMSDterm differences between rigid and flexible dockingterm versus target-cognate side chainterm RMSDsterm, ( b ) average ligand RMSDterm difference for different 1 Å intervals of target-cognate side chainsterm RMSDterm Figure 11 b shows the average RMSDterm difference between rigid and flexible dockingterm for different 1 Å intervals of target-cognate side chainsterm RMSDsterm. For low target-cognate RMSDsterm, which correspond to a highly similar binding siteterm and therefore a situation akin to redockingterm, rigid docking seems to be advantageous on average, as expected. In the other hand, for higher target-cognate side chainsterm RMSDsterm, indicating a lower similarity between binding sitesterm, the situation is less clear. Flexible dockingterm seems to be equivalent or slightly more advantageous, on average, especially at higher target-cognate side chainsterm RMSDsterm. However, as it can be inferred from Fig. 11 a, the number of systems with target-cognate side chainsterm RMSDterm higher than 6 Å is low and therefore the apparent improvement in ligand RMSDterm for flexible dockingterm is inconclusive.

When performing flexible dockingterm, the side chainsterm identified as flexible are included in the calculation of the bounds for the box defining the search space. This results in a larger search space for flexible dockingterm compared to rigid docking, which in turn might require a higher exhaustivenessterm for better samplingterm (although by default the number of Monte Carlo steps is already proportional to the number of degrees of freedom, hence the increased computational cost of flexible dockingterm). It is also worth stressing that we used the default CNN models, which were not explicitly trained for flexible dockingterm.

Given the much higher computational cost of flexible dockingterm, the optimization of Gninaterm default parameters and training of new CNN models for this specific task is outside the scope of the present work and will be addressed in future versions of Gninaterm . However, it is clear that improvements for flexible dockingterm are system-dependent and, therefore, accounting for the increased computational cost, it is reasonable to use rigid docking as the default docking method.

CNN scoringterm performance We evaluate all of the CNN models on a subset of the data that was not seen during training to ensure that the CNN models are able to generalize to unseen protein-ligand systems. We also show the Vinaterm results for the same subset of protein-ligand systems. The CNN models’ docking performance decrease relative to the full sets when looking at the top pose, (comparing Figs. 12 to 4 ). Redockingterm Top1term performance drops from 73% to 68% on the Default Ensembleterm for the filtered set and the full set, respectively, while Vinaterm remains at about 57%. Cross-dockingterm Top1term performance increases from 37% to 42% for the Default Ensembleterm and decreases from 27% to 23% when using the Vinaterm scoring functionterm. Overall, the Default Ensembleterm still shows a strong docking performance boost when used to rescore poses rather than using the Vinaterm scoring functionterm to score poses even on protein-ligand systems not seen during training.

CNN model ensembles evaluated on the subset of proteins and ligands not present in their training datasets. Ensemble models used with the default arguments defined above. TopNterm is the percentage of targets ranked above or at N with a RMSDterm less than 2 Å Finally, we examine the importance of the output CNNscoreterm as a measure of the confidence of the prediction.

We show in Fig. 13 that poses with high CNNscoresterm are more likely to be low RMSDterm to the known binding pose. However, when the CNNscoreterm threshold is close to 1, each CNN ensemble has few poses remaining. Comparing the percentage of complexes remaining when thresholding by CNNscoreterm, we can see that the CNN models are much more confident in the poses when they are performing the redockingterm task: 87% of systems have a top pose with a score higher than 0.8. All of the CNN ensembles identify fewer poses with a high CNNscoreterm for the cross-dockingterm task: only 15% of cross-docked ligands have a pose with a score higher than 0.8 (Additional file 1 : Figure S12). In general, higher CNNscoresterm imply that a pose is more likely to have a low RMSDterm.

Thresholding the top pose by the score determined by the CNN. Top1term is the percentage of top ranked targets with a RMSDterm less than 2 Å Discussion We show that our computational docking software Gninaterm is able to outperform AutoDock Vinaterm scoringterm by using CNN models to rescore generated poses. Without using a CNN model, our software is equivalent to using the Sminaterm docking software, which is a fork of Vinaterm. Gninaterm allows the user to utilize CNN models as scoring functionsterm within the docking pipeline in a variety of ways. The various CNN scoringterm options allow the specified CNN models to replace the scoring functionterm in the samplingterm, refinementterm, and rescoringterm steps of the docking pipeline. We first establish an ensemble of built-in CNN models to be used as the default ensembleterm. This ensemble of models is selected for its docking and runtime performance. The selected ensemble, termed the Default Ensembleterm, is able to exceed the ranking performance of Vinaterm (Top1term increases from 58% to 73% and from 27% to 37% for redockingterm and cross-dockingterm, respectively) while only adding two seconds to the average compute time when utilizing GPUs. It also significantly outperforms any single model or any ensemble of identical models while being significantly faster than an ensemble of all models. The Default Ensembleterm performs almost equally well when performing refinementterm of the sampled poses rather than using the Vinaterm scoring functionterm. Refinementterm with the CNN models is not recommended as the compute cost is significant and the performance is less than simply using the CNN models to rescore the poses refined by the Vinaterm scoring functionterm.

We next derive the default parameters when using the Default Ensembleterm for docking with Gninaterm . autobox_addterm increases the amount of space around the defined binding pocket to allow more volume for samplingterm to investigate. Increases in autobox_addterm decrease the accuracy of the predicted binding pose when the precise binding pocket is known. However, a higher autobox_addterm is necessary if the true binding siteterm is not known. We find that increases in exhaustivenessterm, or the number of Monte Carlo chains run, boost the performance of the docking procedure at the cost of extra computation. We also show that our Default Ensembleterm with an exhaustivenessterm value of 8 (Top1term 73% on redockingterm and 37% on cross-dockingterm) is able to outperform the Vinaterm scoring functionterm with double the exhaustivenessterm (Top1term 58% on redockingterm and 27% on cross-dockingterm). This shows that using the Default Ensembleterm for rescoringterm is more valuable than doubling the amount of samplingterm of the docking procedure. Similarly, increasing the number of poses saved from each Monte Carlo chain, num_mc_savedterm, and the number of poses output by the docking pipeline, num_modesterm, increases the performance of the docking routine at the expense of increased computation. The value of the cnn_rotationterm·var and the value of min_rmsd_filterterm·var do not seem to alter the performance of the docking pipeline. The default arguments for Gninaterm when the binding pocket is explicitly defined are: autobox_addterm (4), exhaustivenessterm (8), num_mc_savedterm (50), num_modesterm (9), min_rmsd_filterterm·var (1.0), and cnn_rotationterm·var (0).

However, when the exact binding pocket is not known and the whole protein is used as the defined binding pocket we see that exhaustivenessterm has a much greater impact on performance. We see that the Default Ensembleterm with exhaustivenessterm 16 (Top1term 43% on redockingterm and 20% on cross-dockingterm) is able to outperform the Vinaterm scoring functionterm even when the exhaustivenessterm is quadrupled to 64 (Top1term 38% on redockingterm and 16% on cross-dockingterm). Again, this shows the importance of using the Default Ensembleterm for rescoringterm rather than increasing exhaustivenessterm with the Vinaterm scoring functionterm. The Default Ensembleterm still shows docking performance improvements with increasing exhaustivenessterm on both redockingterm and cross-dockingterm when the exact binding pocket is unknown. We therefore recommend the largest value of exhaustivenessterm that can be used given time constraints when performing whole protein dockingterm.

We tested the impact of performing cross-dockingterm with flexible side chainsterm (flexible dockingterm). For high binding pocket similarity, flexible dockingterm resulted in a slight deterioration of average ligand top pose RMSDsterm, as expected. However, for low binding pocket similarity there is no clear advantage for flexible dockingterm, whose success is very much system-dependent. In light of this and given the much higher computational cost, rigid docking remains the most suitable default docking procedure.

Finally, we evaluate the ability of the CNN ensembles to score ligand conformationsterm. In order to ensure that our CNN models are generalizing to unseen protein-ligand complexes, we filter the redockingterm and cross-dockingterm benchmark datasets to only include protein-ligand pairs that were not seen during training of the CNNs. We show that our CNN model ensembles are able to outperform the Vinaterm scoring functionterm. The CNN models are able to generalize to unseen examples and properly score the ligand conformationsterm such that the low RMSD posesterm are ranked higher more often than when using the Vinaterm scoring functionterm. The score output can provide a probability of a pose being less than 2 Å from the binding pose. When the CNN outputs a score greater than 0.8 during cross-dockingterm there is at least a 56% probability that the pose is less than 2 Å RMSDterm from the correct pose. If redockingterm is performed, then this probability is 79%. The score output by the CNN models can be used as an indicator of the confidence in the quality of the generated ligand conformationterm; however, particularly when cross-dockingterm, high scoringterm poses are less common (Fig. 13 b).

Gninaterm demonstrates significant improvements over a popular empirical scoring functionterm·def, but many challenges remain. The failure of using gradients from the CNN scoringterm to meaningfully improve docking results points to the need to train new models specifically for this task. Such models, by learning a smoother energy surface than a classical scoring functionterm, may significantly reduce the need for extensive samplingterm and make full CNN docking more computationally tractable. Additionally, new models trained specifically for the flexible dockingterm task would likely improve performance in this area. Gninaterm 1.0 provides a solid framework for these and other improvements. We look forward to receiving contributions and feedback from the community to further enhance the software in future releases as part of our open development processterm.

Supplementary Information Additional file 1: Figure S1 : Analyzing the effect on docking performance when downsampling the cross-docking datasetterm to include a fraction of the protein-ligand pairs per pocket. Figure S2 : Comparison of the RMSDterm(Å) of the top 3 poses output by Gninaterm with no CNN in the docking pipeline and Sminaterm. Figure S3 : Comparison of the minimum and maximum RMSDterm(Å) poses output by Gninaterm with no CNN in the docking pipeline and Sminaterm. Figure S4 : Cross-dockingterm results using a defined binding pocket. Figure S5 : Evaluating different values of cnn emp weight on docking performance when using the Default Ensembleterm. Figure S6 : Time to perform one docking run when using the Default Ensembleterm for “rescore” or refine” in comparison to only using the Vinaterm scoring functionterm. Figure S7 : Evaluating the effect on docking performance when the value of autobox add is changed while using the Default Ensembleterm for rescoringterm. Figure S8 : Evaluating the effect on docking performance when the value of cnn rotation is changed while using the Default Ensembleterm for rescoringterm. Figure S9 : Evaluating the effect on docking performance when the value of min rmsdterm filter is changed while using the Default Ensembleterm for rescoringterm on both the redockingterm and cross-docking datasetsterm. Figure S10 : Cross-dockingterm results using the whole protein as the defined binding pocket. Figure S1 : Ligand RMSDterm distributions for the top pose in the cross-docking datasetterm, for both flexible and rigid docking. Figure S12 : Thresholding the cross-dockingterm results by CNNscoreterm and evaluating Top1term(%) per pocket. A pose is retained if the CNN score is greater than the value on the x-axis. Grey cells indicate that no poses are left in the pocket. Table S1 : Optimal Model Ensemble Selection. Performance given by Top1term, the percent of systems with RMSDterm less than 2 Å from the top pose to the known binding pose. Table S2 : Average time to dock one protein-ligand system from the filtered PDBbind core set v.2016. Comparing runtime when GPU is used to when no GPU is used for docking. Table S3 : Pocket, ligand and receptor identifiers for the systems excluded from the analysis of flexible dockingterm, together with the reason for exclusion.

Acknowledgements The authors thank Monica Dayao, Rosalie Nolen, and Jonathan King for their comments during the preparation of the manuscript. RM thanks Philip Biggin, Andrew Anighoro and Irfan Alibay for fruitful discussions and support. We also thank Zelda, Lily, and Eddie for their canine support.

Author contributions All authors contributed to the development and evaluation of the software. All authors assisted in the preparation of this manuscript. All authors read and approved the final manuscript.

Funding This work is supported by R01GM108340 from the National Institute of General Medical Sciences and is supported in part by the University of Pittsburgh Center for Research Computing through the resources provided. RM is supported by funding from the Biotechnology and Biological Sciences Research Council (BBSRC) National Productivity Investment Fund (NPIF) [BB/S50760X/1] and Evotec (UK) via the Interdisciplinary Biosciences DTP at the University of Oxford [BB/MO11224/1] and was supported by the Google Summer of Code 2019 program.

Availability of data and materials Project name: Gninaterm . Project home page: https://github.com/gninaterm/gninaterm . Operating systems: Linux (Docker container available). Programming language : C++, CUDA. Other requirements: CUDA, Open Babel 3. License: GPL2/Apache License. Any restrictions to use by non-academics: None. All data and scripts used in the creation of this manuscript are available at https://github.com/dkoes/GNINA-1.0 . Source code, scripts, and models are available via the project homepage, https://github.com/gninaterm .All data and scripts used in the creation of this manuscript are available at https://github.com/dkoes/GNINA-1.0 . Source code, scripts, and models are available via the project homepage, https://github.com/gninaterm .

Declarations Competing interests The authors declare that they have no competing interests.

Footnotes References Associated Data Supplementary Materials Additional file 1: Figure S1 : Analyzing the effect on docking performance when downsampling the cross-docking datasetterm to include a fraction of the protein-ligand pairs per pocket. Figure S2 : Comparison of the RMSDterm(Å) of the top 3 poses output by Gninaterm with no CNN in the docking pipeline and Sminaterm. Figure S3 : Comparison of the minimum and maximum RMSDterm(Å) poses output by Gninaterm with no CNN in the docking pipeline and Sminaterm. Figure S4 : Cross-dockingterm results using a defined binding pocket. Figure S5 : Evaluating different values of cnn emp weight on docking performance when using the Default Ensembleterm. Figure S6 : Time to perform one docking run when using the Default Ensembleterm for “rescore” or refine” in comparison to only using the Vinaterm scoring functionterm. Figure S7 : Evaluating the effect on docking performance when the value of autobox add is changed while using the Default Ensembleterm for rescoringterm. Figure S8 : Evaluating the effect on docking performance when the value of cnn rotation is changed while using the Default Ensembleterm for rescoringterm. Figure S9 : Evaluating the effect on docking performance when the value of min rmsdterm filter is changed while using the Default Ensembleterm for rescoringterm on both the redockingterm and cross-docking datasetsterm. Figure S10 : Cross-dockingterm results using the whole protein as the defined binding pocket. Figure S1 : Ligand RMSDterm distributions for the top pose in the cross-docking datasetterm, for both flexible and rigid docking. Figure S12 : Thresholding the cross-dockingterm results by CNNscoreterm and evaluating Top1term(%) per pocket. A pose is retained if the CNN score is greater than the value on the x-axis. Grey cells indicate that no poses are left in the pocket. Table S1 : Optimal Model Ensemble Selection. Performance given by Top1term, the percent of systems with RMSDterm less than 2 Å from the top pose to the known binding pose. Table S2 : Average time to dock one protein-ligand system from the filtered PDBbind core set v.2016. Comparing runtime when GPU is used to when no GPU is used for docking. Table S3 : Pocket, ligand and receptor identifiers for the systems excluded from the analysis of flexible dockingterm, together with the reason for exclusion.

Data Availability Statement Project name: Gninaterm . Project home page: https://github.com/gninaterm/gninaterm . Operating systems: Linux (Docker container available). Programming language : C++, CUDA. Other requirements: CUDA, Open Babel 3. License: GPL2/Apache License. Any restrictions to use by non-academics: None. All data and scripts used in the creation of this manuscript are available at https://github.com/dkoes/GNINA-1.0 . Source code, scripts, and models are available via the project homepage, https://github.com/gninaterm .All data and scripts used in the creation of this manuscript are available at https://github.com/dkoes/GNINA-1.0 . Source code, scripts, and models are available via the project homepage, https://github.com/gninaterm .