Accelerating AutoDock4 with GPUs and Gradient-Based Local Search

Abstract

AutoDock4 is a widely used program for docking small molecules to macromolecular targets. It describes ligand-receptor interactions using a physics-inspired scoring function that has been proven useful in a variety of drug discovery projects. However, compared to more modern and recent software, AutoDock4 has longer execution times, limiting its applicability to large scale dockings. To address this problem, we describe an OpenCL implementation of AutoDock4, called AutoDock-GPU, that leverages the highly parallel architecture of GPU hardware to reduce docking runtime by up to 350-fold with respect to a single-threaded process. Moreover, we introduce the gradient-based local search method ADADELTA, as well as an improved version of the Solis-Wets random optimizer from AutoDock4. These efficient local search algorithms significantly reduce the number of calls to the scoring function that are needed to produce good results. The improvements reported here, both in terms of docking throughput and search efficiency, facilitate the use of the AutoDock4 scoring function in large scale virtual screening.

Fig. 1.

A population processed by a LGA run (Run_ID) can be decomposed into their individuals, and each individual (Ind_ID) can be mapped onto a work-group (WG_ID). The entire set of work-groups is distributed by the GPU runtime scheduler over the available Q compute units (CUs). A CU is a multi-threaded hardware unit capable of processing one work-group (composed of L work-items) at a time. The runs, individuals, and fine-grain tasks are colored according to their associated level of parallelism: high (blue), medium (red), and low (green).

Fig. 2.

The overall AutoDock-GPU workflow consists of a sequence of host (H) and device (D) functions. Program execution always starts and finishes in host functions (depicted at the left side). Time-consuming functions, i.e., kernels, are executed iteratively on the device (depicted at the right side), while their termination is controlled by the host.

Fig. 3.

Distribution of scores returned by LGA runs with increasing number of evaluations (local-search: ADADELTA; local-search rate: 100%). Each violin plot represents scores from 100 LGA runs, and is colored by the input conformation of the ligand used. The global minimum was identified for the protein-ligand complex represented in the upper panel (PDB ID: 1hvy), but not for the complex in the bottom panel (PDB ID: 7cpa) because the distribution of scores did not converge towards a lower bound.

Fig. 4.

Runtime speedups of AutoDock-GPU (on various GPUs) with respect to AutoDock4 on a single CPU core. The number of LGA runs is 100, and the number of evaluations is 2,048,000.

Fig. 5.

Performance scalability (left) and resource-utilization efficiency (right), both in terms of the number of CPU cores. The number of LGA runs is 100, and the number of evaluations is 2,048,000. For each ligand-receptor complex, same program commands were executed on all three CPU instances.

Fig. 6.

Time per evaluation in AutoDock-GPU and AutoDock4. Note that each subplot has a different scale for the y-axis. The red lines are third-degree polynomial fits. These fits were used to interpolate the time per evaluation reported in Table 2.

Fig. 7.

Time per evaluation speedup of AutoDock-GPU (on GTX 1080 Ti and Vega 56 GPUs) over AutoDock4 (on a single CPU core of an AWS c5.18x instance). Note that each subplot has a different scale for the y-axis.

Fig. 8.

Fraction of successful LGA runs as a function of the number of score evaluations, using Solis-Wets and ADADELTA as local-search methods. The upper plots correspond to an easy search problem with only four rotatable bonds (PDB ID: 1tow), while the bottom plots correspond to a moderately difficult ligand with eight rotatable bonds (PDB ID: 1y6b). According to the score criterion, an LGA run is successful if it reports a pose within 1.0 kcal/mol from the global optimum. The black line is a fitted sigmoid curve (Equation 10) that estimates E₅₀, which is the number of evaluations at which 50% of LGA runs are successful. The dashed orange lines are a visual representation of E₅₀ values.

Fig. 9.

Dependency of E₅₀ on the number of ligand rotatable bonds N_rot for Solis-Wets and ADADELTA local-search methods. The LS rate is 6% in the top row, 25% in the middle row, and 100% in the bottom row. The horizontal dashed lines correspond to the lower and upper limits of the number of evaluations used for calculating E₅₀ values.

Fig. 10.

E₅₀ values for ADADELTA and Solis-Wets using 100% LS rate. Each marker represents a protein-ligand complex and is color-coded by the number of rotatable bonds in the ligand. The horizontal and vertical dashed lines correspond to the lower and upper limits of the number of evaluations used for calculating E₅₀ values.

Fig. 11.

Time needed to perform E₅₀ evaluations, using either ADADELTA (x-axis) or Solis-Wets (y-axis). The local search rate was 100%. For a given protein-ligand complex, time-to-E₅₀ is the product of E₅₀ by the corresponding time spent per evaluation. Here, evaluation rates were collected on the Vega 56 GPU platform.