Automatic hyperparameter selection in Autodock

TL;DR

This paper presents an automatic hyperparameter tuning method for AutoDock, using a surrogate-based multi-objective algorithm to improve conformation search efficiency and assist novice users.

Contribution

It introduces a surrogate-based multi-objective optimization approach for automatic hyperparameter tuning in AutoDock, enhancing usability and search performance.

Findings

The method effectively tunes hyperparameters for AutoDock.

Experimental results show improved search efficiency.

The approach is practical and beneficial for users.

Abstract

Autodock is a widely used molecular modeling tool which predicts how small molecules bind to a receptor of known 3D structure. The current version of AutoDock uses meta-heuristic algorithms in combination with local search methods for doing the conformation search. Appropriate settings of hyperparameters in these algorithms are important, particularly for novice users who often find it hard to identify the best configuration. In this work, we design a surrogate based multi-objective algorithm to help such users by automatically tuning hyperparameter settings. The proposed method iteratively uses a radial basis function model and non-dominated sorting to evaluate the sampled configurations during the search phase. Our experimental results using Autodock show that the introduced component is practical and effective.

Figures

Tables

Table 1. TABLE I: Details for the optimized hyperparameters of GA

Hyperparameter	Type	Range	Default
seed	Integer	[0,10000000]	Random
ga_pop_size	Integer	[50,500]	150
ga_elitism	Binary	0,1	1
ga_mutation_rate	Continuous	[0.2,0.99]	0.02
ga_crossover_rate	Continuous	[0.2,0.99]	0.80

Table 2. TABLE II: Details for the optimized hyperparameters of SA

Hyperparameter	Type	Range	Default
seed	Integer	[0,10000000]	Random
tstep	Continuous	[-2.0,2.0]	2.0
qstep	Continuous	[-5.0,5.0]	2.0
dstep	Continuous	[-5.0,5.0]	2.0
rtrf	Continuous	[0.0001,0.99]	0.80
trnrf	Continuous	[0.0001,0.99]	1.0
quarf	Continuous	[0.0001,0.99]	1.0
dihrf	Continuous	[0.0001,0.99]	1.0
accs	Integer	[100,30000]	30000
rejs	Integer	[100,30000]	30000
linear_schedule	Binary	0,1	1

Table 3. TABLE III: Details for the optimized hyperparameters of the adopted local search

Hyperparameter	Type	Range	Default
seed	Integer	[0,10000000]	Random
sw_max_its	Integer	[100,1000]	300
sw_max_succ	Integer	[2,10]	4
sw_max_fail	Integer	[2,10]	4

Table 4. TABLE IV: The average docking results of the SA, GA, LS and HB algorithms after 10 runs for the first case study.

Algorithm	Energy	RMSD
Tuned ${\text{SA}}_1^{\ast }$	-13.4990	4.0450
Tuned ${\text{SA}}_2^{\ast }$	-10.9280	3.4140
${\text{SA}}_1$	-10.9260	3.5460
${\text{SA}}_2$	-10.9260	3.5460
Tuned ${\text{GA}}_1^{\ast }$	-13.7110	3.4960
Tuned ${\text{GA}}_2^{\ast }$	-11.9070	2.4550
${\text{GA}}_1$	-10.3610	5.2730
${\text{GA}}_2$	-10.1120	3.2880
Tuned ${\text{LS}}_1^{\ast }$	-17.4270	0.2190
Tuned ${\text{LS}}_2^{\ast }$	-17.1530	0.0860
${\text{LS}}_1$	-17.2990	0.1980
${\text{LS}}_2$	-17.1920	0.1400
Tuned ${\text{HB}}_1^{\ast }$	-13.7110	3.4960
Tuned ${\text{HB}}_2^{\ast }$	-11.9070	2.4550
${\text{HB}}_1$	-12.6610	2.6910
${\text{HB}}_2$	-12.6610	2.6910

Table 5. TABLE V: The average docking results of the SA, GA, LS and HB algorithms after 10 runs for the second case study.

Algorithm	Energy	RMSD
Tuned ${\text{SA}}_1^{\ast }$	-15.0124	6.2430
Tuned ${\text{SA}}_2^{\ast }$	-13.0790	6.0120
${\text{SA}}_1$	-14.0790	6.5380
${\text{SA}}_2$	-12.9610	6.4640
Tuned ${\text{GA}}_1^{\ast }$	-12.728	6.722
Tuned ${\text{GA}}_2^{\ast }$	-12.728	6.722
${\text{GA}}_1$	-14.0030	6.6350
${\text{GA}}_2$	-12.1460	6.4200
Tuned ${\text{LS}}_1^{\ast }$	-15.6910	6.6340
Tuned ${\text{LS}}_2^{\ast }$	-13.4850	6.4690
${\text{LS}}_1$	-14.8500	6.6480
${\text{LS}}_2$	0.1520	6.3570
Tuned ${\text{HB}}_1^{\ast }$	-16.9080	6.6600
Tuned ${\text{HB}}_2^{\ast }$	-15.6030	6.4730
${\text{HB}}_1$	-16.0210	6.7040
${\text{HB}}_2$	-15.8590	6.6160

Equations

$$\mathbf{S}={\left[{\mathbf{x}}^{(1)},...,{\mathbf{x}}^{(n)}\right]}^{\text{T}}\in {\mathbb{R}}^{n\times m},\mathbf{x}=\left\{x_1,...,x_m\right\}\in {\mathbb{R}}^m$$

$$\mathbf{y}={\left[y^{(1)},...,y^{(n)}\right]}^{\text{T}}={\left[y({\mathbf{x}}^{(1)}),...,y({\mathbf{x}}^{(n)})\right]}^{\text{T}}\in {\mathbb{R}}^n$$

$$E_{inter}=Q_{bound}^{R-L}+Q_{unbound}^{R-L}$$

$$\begin{gathered} Q=W_{vdw}\sum _{i,j}\left(\frac{A_{ij}}{r_{ij}^{12}}-\frac{B_{ij}}{r_{ij}^6}\right)+ \\ {W}_{hbond}\sum _{i,j}E(t)\left(\frac{C_{ij}}{r_{ij}^{12}}-\frac{D_{ij}}{r_{ij}^{10}}\right)+ \\ {W}_{elec}\sum _{i,j}\frac{q_i{q}_j}{\epsilon (r_{ij})r_{ij})}+W_{sol}\sum _{i,j}(S_j{V}_i+S_i{V}_j)e^{(-\frac{r_{ij}^2}{2{\sigma }^2})} \end{gathered}$$

$$\tilde{f}(\hat{\mathbf{x}})=f(\hat{\mathbf{x}})+ϵ$$

$$\tilde{f}(\hat{\mathbf{x}})=\sum _{i=1}^n{w}_{i}h(\left|\hat{\mathbf{x}}-{\mathbf{x}}_i\right|)+{\mathbf{b}}^{\top }\hat{\mathbf{x}}+a$$

$$\left(\begin{array}{cc}\mathbf{\Phi }& \mathbf{P}\\ {\mathbf{P}}^{\top }& 0\end{array}\right)\left(\begin{array}{c}\mathbf{w}\\ \mathbf{c}\end{array}\right)=\left(\begin{array}{c}{f}_{\psi }\\ 0\end{array}\right)$$

$${\mathbf{P}}^{\top }=\left(\begin{array}{cccc}{\mathbf{x}}_1& {\mathbf{x}}_2& \cdots & {\mathbf{x}}_n\\ 1& 1& \cdots & 1\end{array}\right)$$

$${\mathbf{\sigma }}_j={\mathbf{x}}_{best}+\gamma \otimes \text{v},j=\mathrm{1...}N$$

Full text

I Introduction

AutoDock is a protein-ligand docking tool which has now been distributed to more than 29000 users around the world[1]. The application for this software arises from problems in computer-aided drug design. The ideal procedure would optimize the interaction energy between the substrate and the target protein. It employs a set of meta-heuristic and local search optimization methods to meet this demand. The Autodock puts protein-ligand docking at the disposal of non-expert users how often do not know how to choose among the dozens of configurations for such algorithms. This motivated us to bring the benefits of hyperparameter tuning to users of AutoDock.

The global optimization problem arises in many real-world applications as well as molecular docking. A standard continuous optimization problem seeks a parameter vector $\mathbf{x}^{*}$ to minimize an objective function $f(\mathbf{x}):\mathbb{R}^{D}\rightarrow\mathbb{R}$, i.e. $f(\mathbf{x}^{*})<f(\mathbf{x})$ for all $\mathbf{x}\in\Omega$, where $\Omega=\mathbb{R}^{D}$ denotes the search domain (a maximization problem can be obtained by negating the $f)$. Over the years, optimization algorithms have been effectively applied to tackle this problem. These algorithms, however, may need to be fine-tuned which is a major challenge to their successful application. This is primarily due to the highly nonlinear and complex properties associated with the objective function. In this context, the advantages of algorithm configuration techniques become clear.

The hyperparameter tuning domain has been dominated by model-based optimization [2] that adopt probabilistic surrogate models to replace in part the original computationally expensive evaluations. They construct computationally cheap to evaluate surrogate models in order to provide a fast approximation of the expensive fitness evaluations during the search process. Surrogate based techniques try to model conditional probability $p(y|\varphi)$ of a $m-$dimensional configuration $\varphi$ given $n$ observations S:

[TABLE]

with the corresponding evaluation metrics $\mathbf{y}$:

[TABLE]

The essential question arises in model based algorithms is which individuals should be chosen to be evaluated using the exact fitness function. This is most characterized by making a good balance between the exploration and exploitation capabilities. One of the earliest researches in this direction was performed by Jones et al. who adopted a Kriging surrogate to fit the data and make a balance between exploration and exploitation [3]. The exploration property of their proposed efficient global optimization (EGO) algorithm is enhanced by the fact that the expected improvement (EI) (an acquisition function) is conditioned on points with large uncertainty and low values of the surrogate. In another study, Booker et al. [4] sought out a balanced search strategy by taking into account sample points with low surrogate predictions and high mean square error. Moreover, Wang et al. [5] introduced the mode-pursuing approach which favors the trial points with low surrogate values according to a probability function. Regis and Shoemaker [6] also put forward an approach according to which the next candidate point is chosen to be the one that minimizes the surrogate value subject to distance from previously evaluated points. The distance starts from a high value (global search) and ends with a low value (local search). They also proposed Stochastic Response Surface (SRS) [7] algorithm that cycles from emphasis on the objective to emphasis on the distance using a weighting strategy. The SRS, however, mitigated some of the requirements for inner acquisition function optimization. For this reason, we focused on proposing a new algorithm configuration approach based on SRS model. The proposed algorithm, called as MO-SRS, modify the SRS so as to be able to handle multi-objective problems. More precisely, MO-SRS is equipped with the idea of multi-objective particle swarm optimization (PSO)[8]. We used the MO-SRS to make a balance between intermolecular energy and the Root Mean Square Deviation (RMSD) during the hyperparameter optimization in Autodock.

The rest of the paper is organized as follows. Section II provides us with a brief review on the related works. Section III gives a brief description for the docking problem. Section IV elaborates technical details of our proposed approach. In Section V, the performance of the introduced components is investigated by conducting a set of experiments. The last section summarizes the paper and draws conclusions.

II Related works

Sequential Model-based Algorithm Configuration (SMAC) [2], Hyperband [9], Spearmint [10], F-race [11] and Tree-structure Parzen Estimator (TPE) [12] are examples of well known methods for hyperparameter optimization. SMAC adopted a random forests model and Expected Improvement (EI) to compute $p(y|\varphi)$. Similarly, TPE et al. defined a configuration algorithm based on tree-structure Parzen estimator and EI. To tackle the curse of dimensionality, TPE assigns particular values of other elements to the hyperparameters which are known to be irrelevant. Ilievski et al. [13] proposed a deterministic method which employs dynamic coordinate search and radial basis functions (RBFs) to find most promising configurations. By using the RBFs [14] as surrogate model, they mitigated some of the requirements for inner acquisition function optimization. In another work [15], the authors put forward neural networks as an alternative to Gaussian process for modeling distributions over functions. Interestingly, Google introduced Google Vizier [16], an internal service for surrogate based optimization which incorporates Batched Gaussian Process Bandits along with the EI acquisition function.

Although the above-mentioned approaches have been proven successful, they are not able to handle several objective functions during the tuning process. Consequently, this paper presents a novel multi-objective algorithm which integrates the idea of a multi-objective meta-heuristics with the SRS method to find a subset of promising configurations.

III Molecular docking

Molecular docking enables us to find an optimal conforma- tion between a ligand $L$ and receptor $R$. In another word, it predicts the preferred orientation of $L$ to R when bound to each other in order to form a stable complex. A schematic example of this procedure is illustrated in Fig.1. This problem can be formulated as a multi-objective problem consisting of minimizing the RMSD score and the binding energy $E_{inter}$. In Autodock, energy function $E_{inter}$ is defined as follows [1]:

[TABLE]

[TABLE]

In (3), $Q_{bound}^{R-L}$ and $Q_{unbound}^{R-L}$ represent the states of ligand-protein complex in bound and unbound modes, respectively.

The pairwise energetic terms take into the account evaluations for repulsion ($vdw$), electrostatics ($elec$), desolvation ($sol$) and hydrogen bonding ($hbond$). The constant weights $W_{vdw}$, $W_{elec}$, $W_{hbond}$ and $W_{sol}$ belong to Van der Waals, electrostatic interactions, hydrogen bonds and desolvation, respectively. Moreover, $r_{ij}$, $A_{ij}$, $B_{ij}$, $C_{ij}$ and $D_{ij}$ denotes Lennard-Jones parameters, while $E(t)$ function is the angle-dependent directionality. Also, $V$ is the volume of atoms that surround a given atom, weighted by a solvation parameter ($S$). For a detailed report of all the variables please refer to [17].

IV The proposed method

This section discusses in detail the main components of the proposed MO-SRS. We extend the idea of stochastic RBF to the multi-objective case to be suitable for the Autodock application. Compared to the evolutionary algorithms like genetic algorithm, MO-SRS need less computational time by virtue of surrogate modeling techniques. A workflow of the introduced approach is illustrated in Fig 2.

The MO-SRS starts the optimization procedure by generating a set of random configurations (i.e., initial population). During this phase, it might be possible to miss a considerable portion of the promising area due to the high dimensionality of the configuration space. It should be noticed that we have a small and a fix computational budget and increasing size of the initial population cannot remedy the issue. Furthermore, it is crucial for a model-based algorithm to efficiently explore the search space so as to approximate the nonlinear behavior of the objective function. The Design of Computer Experiments (DoCE) methods are often used to partially mitigate high dimensionality of the search space. Among them, MO-SRS utilized the Latin Hypercube Sampling (LHS) [18] which is a DoCE method for providing a uniform cover on the search space using a minimum number of population. The main advantage of LHS is that it does not require an increased initial population size for more dimensions.

As the next step, we evaluate all the generated configurations $\mathbf{x}_{i}(i=1,2,\cdots,n)$ using the Autodock to yield outcomes $y_{i}^{(1)}=f^{(1)}(\mathbf{x_{i}})$ and $y_{i}^{(2)}=f^{(2)}(\mathbf{x_{i}})$. Here, ${f_{1}}$ and ${f_{2}}$ denote the energy and RMSD function variables, respectively. Thereafter, we adopt surrogate model techniques to approximate the fitness function using this data set.

Surrogate models offer a set of mathematical tools for predicting the output of an expensive objective function $f$. Particularly, they try to predict the fitness function $f$ for any unseen input vector $\hat{x}$ according to computed data points $(x_{i},y_{i})$. Given solution $\boldsymbol{\hat{x}}$ and objective function $f$, a surrogate model $\hat{f}$ can be defined as in (4), where $\epsilon$ is the approximation error.

[TABLE]

Among different surrogate models, we utilized the Radial Basis Function (RBF) which is a good model for high dimensional problems [19, 14]. It belongs to an interpolation method for scattered multivariate data which considers all the sample points. In this light, it introduces linear combinations based on a RBF $h(\mathbf{x})$ to approximate the desired response function $f$; as presented in (5).

[TABLE]

In (5), $w_{i}$ shows $i^{th}$ unknown weight coefficient, $h$ determines a radial basis function, $\hat{\mathbf{x}}$ is an unseen point and $\epsilon$ are independent errors with a variance $\sigma^{2}$. A radial function $h:\mathbb{R}^{m}\rightarrow\mathbb{R}$ has the property $h(\mathbf{x})=h(\left\|x\right\|)$. By having a suitable kernel $h$; the unknown parameters $a,\mathbf{b}$ and $\mathbf{w}$ could be obtained by solving the following system of linear equations:

[TABLE]

Here, $\mathbf{\Phi}$ is a $n\times n$ matrix with $\Phi_{i,j}=h(\left\|\mathbf{x}_{i}-\boldsymbol{x}_{i}\right\|)$, $\mathbf{c}=(\mathbf{b}^{\top},a)^{\top}$ and

[TABLE]

The linear system (6) has a unique solution and it can be used to approximate function values.

As opposite to the standard SRS method, we have two main objective functions which should be optimized: the intermolecular energy and the RMSD. Hence, we repeat the above mentioned procedure for both the objective functions. More precisely, we will train two surrogate models $\hat{f_{1}}$ and $\hat{f_{2}}$ to approximate the energy and RMSD output variables, respectively.

Now we are ready to use the trained surrogate models inside a pre-selection strategy. Accordingly, we generate a set of $N$ candidates to be evaluated using $\hat{f_{1}}$ and $\hat{f_{2}}$. In the single-objective SRS code, the following equation generates neighborhoods using the current best solution $\hat{f_{2}}$, a randomly produced vector $\textup{{v}}\in\mathbb{R}^{m}$ and adaptive parameter $\gamma$:

[TABLE]

Generally speaking, in multi-objective optimization the best solutions need to be chosen based on two objectives and there is no single best solution. So, we cannot directly apply the standard search operators of the SRS. Moreover, we have to update our surrogate models $\hat{f_{1}}$ and $\hat{f_{2}}$ based on a best obtained solution which is another problem for the MO-SRS. Having this in mind, we borrowed the idea of leader selection in multi-objective PSO [8] algorithm which handles the same problem (update equation of the PSO also depends on the global best solution and finding a leader particle is an important task). Inferred from [8], MO-SRS first builds an archive which contains the Pareto-optimal solutions. A solution vector is Pareto-optimal if there is not another solution that dominates it. In mathematical terms, solution $x_{1}$ is said to dominate another solution $x_{2}$, if:

• •

$f_{i}(x_{1})\leqslant f_{i}(x_{2})$ for all $i\in\left\{1,2,...,k\right\}$ and

• •

$f_{i}(x_{1})\leqslant f_{i}(x_{2})$ for at least one $i\in\left\{1,2,...,k\right\}$

The MO-SRS divides the objective space into hypercubes and assigns a fitness value to each hypercube based on the number of Pareto-optimal solutions that lie in it. Then, it employs the roulette-wheel to find the superior hypercube. A randomly selected solution from that hypercube is determined to be the best solution $\boldsymbol{x}_{best}$. All the $N$ candidates are evaluated by the previously trained surrogates $\hat{f_{1}}$ and $\hat{f_{2}}$. Again, we select the best solution according to the above mentioned strategy. This best solution is then used to update our surrogate models. All the aforementioned steps continue until some stopping criteria are met. To sum up, MO-SRS approach is illustrated in Fig. 3. The MO-SRS will be available soon on https://github.com/ML-MHs/Auto-Autodock.

Fig. 3. Pseudocode of the MO-SRS method

V Experimental results

In this section, the performance of MO-SRS hyperparameter tuning method on Autodock task is investigated. We considered two scenarios from Autodock documentations which should be minimized: 1)1DWD and 2) HSG1. All simulations are performed using Matlab. We adopted the same implementation and parameter configuration for the SRS as suggested in [20]. The Cubic function is used as the kernel in the RBF.

We applied the MO-SRS to tune the hyperparameters of Autodock for the genetic algorithm (GA) [21], simulated annealing (SA) [22], a local search (LS) called as pseudo-Solis-Wet [1] and a hybrid method (HB) which combines the GA and LS. The GA is a well-known evolutionary optimization algorithm which is highly sensitive ot the initial value of its parameters. The performance of the GA is dependent on the structure of the problem at hand and do not scale well with complexity. Consequently, parameters of GA like crossover, population size and mutation should be chosen with care. A small population size will lead to the early convergence problem, while using a large population increases the computational cost. The configuration tuning of GA becomes even more challenging when there is a correlation between its hyperparameters. For example, mutation is more effective on smaller population sizes, while the crossover is likely to benefit from large populations. All the mentioned reasons make the GA a suitable algorithm for benchmarking the performance of MO-SRS. The same situation could happen for the considered LS and SA algorithms. From other points of view, the adopted algorithms are introduced to measure search performance of the MO-SRS under different dimensions. The LS and GA belong to low dimensional problems, SA is a medium dimension and HB algorithm is a high dimensional hyperparameter tuning problem. A detailed information for the optimized configurations for GA, SA and LS are presented in Tables I-III, respectively.

The obtained results are reported in Tables IV and V. The MO-SRS should obtain a reliable search performance using a limited computational budget and so the number of evaluations is set to 100. For each algorithm, the search bounds are reported in Tables I-III. To reduce the influence of stochastic error, experiments are repeated 10 times for each problem. In Tables IV-V, the *tuned * prefix denotes the optimized algorithms using the MO-SRS. The MO-SRS offers a set of final solutions for each of the cases and we used a Pareto-optimal graph to illustrate such solutions; as depicted in Fig. 4. However, it should be noticed that the results offered by the default methods in Autodock are based on single-objective algorithms. For this reason, we compared the obtained results of the MO-SRS according to each of those objectives. The subindex 1 in tables IV-V shows the performance of the compared algorithms for the binding energy and subindex 2 denotes the same for the RMSD criterion. The results are averaged over 10 runs. The energy unit is in kcal/mol and RMSD unit is in Angstroms.

The results of the first case study are presented in Table IV. As it can be seen, MO-SRS converges closer to the global optimum. In the case of LS method, we can see that the tuned algorithm yields a minimized energy -17.1530 and RMSD 0.0860. In docking domain, RMSD $<$ 1.5 Angstrom is always wise to consider. Similarly, MO-SRS provides more accurate results for the second scenario according to the results in Table V.

In these tables, we can see how the introduced method tries to find more than one objective simultaneously. In this regard, there are two main points. The first one is that the diversity of the obtained solutions strongly depends on the performance of the algorithm at hand. For example, LS performs better than the GA and consequently we can see the Pareto-optimal for LS contains more solution. The second point is that all Pareto-optimal solutions are obtained only after 100 evaluations. This confirms our hypothesis that MO-SRS can effectively be use to tune hyperparameters of the Autodock. The main components of MO-SRS are the adopted multi-objective and surrogate modeling techniques which can be further improved in the future works.

VI Conclusion

A typical docking scenario in Autodock includes applying optimization algorithms. Following those steps, users should select a set of appropriate hyperparameters to maximize the performance of their final results. As it is often beyond the abilities of novice users, we proposed to develop a multi-objective approach for automatically tuning highly parametric algorithms of the Autodock. Automating end-to-end process of applying the proposed MO-SRS offers the advantages of more accurate solutions. The experimental results clearly show that the obtained results by MO-SRS outperform hand-tuned models. We hope that the introduced MO-SRS helps non-expert users to more effectively apply Autodock to their applications.