=Paper=
{{Paper
|id=None
|storemode=property
|title=A Comparative Study of Different Optimization Algorithms for Molecular Docking
|pdfUrl=https://ceur-ws.org/Vol-819/paper5.pdf
|volume=Vol-819
|dblpUrl=https://dblp.org/rec/conf/iwsg/AfanasievOPRSS11
}}
==A Comparative Study of Different Optimization Algorithms for Molecular Docking==
https://ceur-ws.org/Vol-819/paper5.pdf
3rd International Workshop on Science Gateways for Life Sciences (IWSG 2011), 8-10 JUNE 2011
A Comparative Study of Different Optimization Algorithms for Mo-
lecular Docking
Alexander Afanasiev1 , Igor Oferkin3, Mikhail Posypkin1,* , Anton Rubtsov1, Alexey Sulimov3
and Vladimir Sulimov2,3
1
Institute for Systems Analysis RAS, Address 117312, Moscow, pr. 60-letiya Oktyabrya, 9, mposypkin@mail.ru
2
Research Computer Center at Lomonosov Moscow State University, 119992, Leninskie Gory 1, bldg 4, Moscow 11999 2,
Moscow, Russia
3
Dimonta, Ltd, 117186, Nagornaya Str. 15, build. 8, Moscow, Russia
ABSTRACT Having fast and robust optimization algorithms for solving the
Motivation: Computer modeling of protein-ligand interactions is one problem (1) is crucial for an efficient docking. In this paper we
of the most important phases in a drug design process. The core evaluated different optimization algorithms for resolution of the
part of this modeling is a resolution of a global unconstrained optimi- problem (1). We performed numerical experiments, analyzed the
zation problem. This paper presents a comparative computational results and suggested the most successful combination of optimiza-
experiments aimed at studying the efficiency of the different optimi- tion algorithms for this problem.
zation methods applied to the docking problem. We present experi-
mental results for different optimization algorithms and draw conclu-
sions about their efficiency.
2 RELATED WORK
The comparative efficiency of different optimization algorithms
1 INTRODUCTION has been studied in various papers (Rosin 1997, Morris 1998, Ta-
Many diseases are caused by foreign protein activity or protein vares 2008, Tavares 2009). Authors consider global optimization
malfunction. For the treatment of these diseases we can try to approaches (genetic algorithms, simulated annealing), local search
block these proteins by small organic molecules. These molecules techniques (L-BFGS, Solis-Wets method) and their combination.
can selectively bind to proteins and thus block their work. This The best results were obtained by a combined approach when a
simplified conception allows development of the drugs, using pur- genetic algorithm is combined with a local search. In this approach
poseful design of new organic compounds – inhibitors for the giv- a fraction of all individuals in a generation are further optimized by
en target-protein. The selection of the effective ligands for inhibi- applying a local search technique.
tion of the target enzyme is usually very laborious, long, and ex-
pensive process. Contemporary molecular modeling tools can ac- Two local search techniques addressed in the literature showed the
celerate this process and make it much less expensive. Virtual best results. The first method proposed in (Solis 1981) is a direct
screening by means of ligands docking is widely recognized as a search method with an adaptive step size, which performs a ran-
helpful approach in modern drug design (Kitchen 2004, Zoete domized local minimization of a given candidate solution. Depend-
2009). ing on whether a new solution if found or not, a success or a failure
is recorded. If several successes occur in a row, the step is in-
The goal of docking is to find the positions of interacting ligand creased to move more quickly. If the opposite occurs, the step is
and a protein with a minimal energy. The stability of this position decreased. A bias term is applied to drive the search in successful
is characterized by the energy value. The ligands with minimal directions. The method terminates when a certain lower-bound
values are candidates for further consideration as potential drugs. threshold for is passed or when a maximum number of steps is
Thus the docking problem is reduced to the global optimization reached. The second method is the Broyden-Fletcher-Goldfarb-
problem Shanno method (Nocedal 2006). L-BFGS is a quasi-Newton me-
thod, where both the function to minimize and its gradient must be
→ , ∈ (1) supplied by the user. The method stops as soon as it finds a local
optimum or when a threshold number of iterations is exceeded.
where is a tuple that determines the ligand position and is a
bounding box restricting the search area to a reasonable region. The direct evaluation of the protein-ligand interaction energy is
computationally expensive. Therefore in modern docking tools the
grid of potentials representing protein-ligand interactions is calcu-
lated separately before the docking procedure and the energy is
approximated using precalculated values. The resulting function is
*
To whom correspondence should be addressed.
Copyright © 2011 for the individual papers by the papers' authors. Copying permitted only for private and academic purposes. This volume is published and copyrighted by its editors.
3rd International Workshop on Science Gateways for Life Sciences (IWSG 2011), 8-10 JUNE 2011
not differentiable in the grid nodes and thus L-BFGS method is not Bond lengths and valence angles have been frozen in the
applicable. course of the docking procedure.
However of the energy is calculated directly local search methods
employing gradient information can be advantageous. It is shown 4 OPTIMIZATION METHODS UNDER TEST
in (Tavares 2009) that L-BFGS method (Nocedal 2006) signifi- We considered three types of optimization methods: local, semi-
cantly improves the performance of genetic algorithms and gives local and global.
better results w.r.t. Solis-Wets method.
4.1 Local methods
Though the efficiency of various search methods was ad-
dressed in the literature some important methods were The goal of local methods is to find the local minimum i.e. the
point that gives the minimal function value in a some neighbor-
missed. For instance semi-local methods e.g. Monotonic
hood. We considered four methods described in table 1.
Basin Hopping that proved to be very efficient for atomic
cluster conformation problem (Leary 2000) were not consi-
dered at all. In this paper we classify search methods into Table 1. Local optimization methods under test
three groups: local, semi-local and global and performs the
systematic evaluation of the several techniques and their Name Reference Brief description
combinations.
CG (Polak 1971) Conjugate Gradients Method
TNC (Nash 2000) Truncated Newton Method
3 CALCULATING PROTEIN-LIGAND LBFGS (Nocedal 2006) The limited memory Broyden–
INTERACTION ENERGY Fletcher–Goldfarb–Shanno
(BFGS) method
Powell (Powel 1964) Powell's method
The protein-ligand interaction energy is the objective function in
the problem (1). Its accurate evaluation is crucial for a success of
the whole interaction simulation and thus for the validity of the The conjugate gradient method is a seminal optimization method
numerical experiments. that is explained in almost every global optimization textbook. It
uses conjugate directions instead of the local gradient for going
Two very popular programs implementing docking algorithms downhill. We used the Polak-Ribere form for calculating conjugate
GOLD (Cole 2005) and AutoDock (Morris 1998, Morris 2005) directions.
employ too simplified force field, either neglecting electrostatic
interaction in GOLD, or too simplified treating of desolvation The classical Newton method requires fewer iterations than conju-
terms in AutoDock. gate gradients but each iteration involves the resolution of the sys-
tem of linear algebraic equations. The truncated Newton methods
For our study we considered the interaction model called SOL are based on the idea that an exact solution of the Newton equation
proposed in (Romanov 2004, Romanov 2008, Oferkin 2011). The at every step is unnecessary and can be computationally wasteful
main idea of this model is to describe with maximal possible accu- in the framework of a basic descent method. Any descent direction
racy the protein-ligand interactions, using the docking procedure will suffice when the objective function is not well approximated
based on contemporary molecular mechanics. The main distinctive by a convex quadratic and, as a solution to the minimization prob-
features of SOL are: lem is approached, more effort in solution of the Newton equation
• A rigid target-protein with the active site represented by may be warranted.
a set of grids for different type potentials, describing pro-
tein-ligand interactions (electrostatic, Van der Waals Another very successful approach to local unconstrained optimiza-
(VdW) forces) in the frame of Merck Molecular Force tion is Broyden–Fletcher–Goldfarb–Shanno (BFGS) method. This
Field (MMFF) (Halgren 1996). algorithm belongs to the family of quasi-Newton methods. In these
• Quite rigorous description of solvation-desolvation ef- methods the Hessian matrix of second derivatives need not be eva-
fects upon ligand binding process, based on Generalized luated directly. Instead, the Hessian matrix is approximated using
Born approximation (Ghosh 1998) and included in the rank-one updates specified by gradient evaluations (or approximate
set of potential grids. gradient evaluations). The L-BFGS algorithm studied in this paper
• The grid of potentials representing protein-ligand inte- stores only a few vectors that represent the approximation implicit-
ractions are calculated separately before the docking pro- ly. Due to its moderate memory requirement, L-BFGS method is
cedure. Electrostatic, VdW and solvation-desolvation po- particularly well suited for optimization problems with a large
tentials were calculated on the 101x101x101 grid inside number of variables.
this cube
• All ligands are considered as fully flexible, i.e. all topo- The three methods outlined above require the evaluation of the
logically available torsion degrees of freedom were un- objective function gradient. However for docking problem the
frozen and allowed to rotate freely, directed only by li- gradient is either not known symbolically or computationally in-
gand internal energy preferences in the frame of MMFF.
3rd International Workshop on Science Gateways for Life Sciences (IWSG 2011), 8-10 JUNE 2011
tractable. Thus one has to approximate the gradient numerically or sequence of uniformly distributed random generated points in a
use methods that don’t rely on gradient information. For the Pow- feasible box. After that each point is used a starting point for semi-
ell’s method the objective function need not be differentiable, and local methods.
no derivatives are taken. The method minimizes the function by a
bi-directional search along each search vector. The new position
can then be expressed as a linear combination of the search vec- 5 EXPERIMENTS
tors. The new displacement vector becomes a new search vector, We performed several experiments for different protein-ligand
and is added to the end of the search vector list. Meanwhile the pairs. Methods demonstrate the same relative behavior for all va-
search vector which contributed most to the new direction, i.e. the riants and thus we present data only for one pair: the target protein
one which was most successful, is deleted from the search vector is thrombin (PDB code 1O2G) and the ligand is 4-aminopyridine
list. The algorithm iterates an arbitrary number of times until no in the protonated form.
significant improvement is made. The method is useful for calcu-
lating the local minimum of a continuous but complex function,
especially one without an underlying mathematical definition, 5.1 Testing results for local methods
because it is not necessary to take derivatives.
The results for local methods are summarized in the table 3. The
average energy is calculated as a mean value of 64 runs with ran-
4.2 Semi-local methods domly generated initial points. The best energy is a lowest value
For functions with multiple extremes a local search methods get found throughout these runs. The percentage of errors shows the
stuck in local minima. The semi-local methods can “escape” from number of runs that produced a value greater than the initial value.
a local minimum by exploring its neighborhood. Such methods
have proved their efficiency for molecular conformation problems
Table 3. Testing results for local methods
(Wales 1997). We considered two such methods summarized in the
Table 2.
Name Avrg. Energy Best Energy Errors (%) Avrg. Time(s)
Table 2. Semi-Local optimization methods under test
Name Reference Brief description Initial 15385,25 4728,66 0,00 0,00
CG 5790,28 874,80 9,52 1,81
MBH (Leary 2000) Monotonic Sequence Basin Hopping TNC 4222,24 1322,46 0,00 3,90
BP (Panteleev 2005) Best Probe Method LBFGS 3119,95 782,21 1,59 5,92
Powell 611,49 28,66 0,00 5,12
The monotonic sequence basin hopping method tries to improve
the minimum until the number of attempts exceeds the threshold As expected the local methods generally improve the energy value.
value . Starting from some point ∈ MBH performs the fol- The Powell method remarkably outperforms methods that rely on
lowing steps ( = 0 at the beginning): Taylor formula. This is an expectable result as the objective func-
tion obtained as a result of piece-wise linear approximation on a
1. Select randomly a point ∈ . mesh is non-differentiable in minima. Therefore the Taylor series
2. Obtain point by applying local minimization to gives a poor approximation for a goal function in a neighborhood
∈ : = . of such points. Thus the gradient information can only be used to
3. If ≤ then assign ≔ , ≔ 0. Otherwise as- bias the search to the descending direction but not as a stopping
sign ≔ + 1.. criterion.
4. If ≥ . then finish, otherwise go to 1.
where = ∈ : | − | ≤ !, = 1, … , # is a box
5.2 Testing results for semi-local methods
neighborhood of the given radius ! .
Thus the MBH method is parameterized by the threshold value , Results of the previous section clearly indicate that the Powell
the neighborhood radius ! and the local search method . method is the best local search strategy among the considered set.
Thus we used this method as a local method in MBH and BP
The best probe method is based on the same idea as MBH. But semi-local methods. After a set of experiments we found that the
unlike MBH it chooses points on the n-dimensional sphere rather best results are obtained if MBH uses the radius ! = 0.1 and the
than in a box and uses adaptive neighborhood size. Starting with a BP uses the initial sphere radius ! = 0.5 . The accuracy of results
large radius it reduces its size if attempts didn’t lead to an im- for both methods depends on the threshold values: the higher thre-
provement. shold value the lower the minimum. To put both methods in the
equal conditions we set the threshold value to 30 and 90 for BP
and MBH respectively. With such parameters both methods take
4.3 Global methods approximately the same time.
The goal of global methods is to diversify the search in order to The Table 4 compares the results of the basin hopping and best
explore the whole feasible set. At the moment we considered only probe methods. The results were averaged over 10 runs with dif-
one globalization strategy: the Monte-Carlo method. It generated a ferent random initial points generated in a box .
3rd International Workshop on Science Gateways for Life Sciences (IWSG 2011), 8-10 JUNE 2011
Table 4. Testing results for semi-local methods optimization problems in the grid and supercomputers. This tool is
currently ported to the desktop grid and is validated for running in
the EDGeS VO. Docking will be one of the supported applications.
Name Avrg. Energy Best Energy Avrg. Time(s)
Another one – atomic cluster conformation problem (Leary 2000)
is already running in production. With such approach we’ll reuse
Initial 16885,3 5628,76 0,00
the one deployed application for different problems. This is very
MBH -73,44 -134,60 451,37 important for desktop grids and combined infrastructures where the
BP -27,77 -131,93 491,84 deployment and validation requires lots of efforts.
Both methods under test gives approximately the same best results, The docking is only one element of a complex workflow in a drug
however the average behavior of MBH is much better. design process. The good software environment is crucial for mak-
ing the docking software tools useful for a wide community of
researchers. One of such environments is described in (Kim 2008).
5.3 Testing results for global methods It flexibly integrates the convenient Web-based user interface and
At the moment we tried only one globalization strategy: Monte the powerful processing back-end deployed in the Grid. Such soft-
Carlo method that generates random initial points for MBH me- ware architecture seems to be the standard approach for docking
thods. Points are generated in the bounding box . Table 5 summa- and other applications demanding for the huge computing power.
rizes results obtained from 10 runs of this combination. The thre- To implement a consistent and flexible software environment for
shold value for MBH method was set to 90 and the initial radius drug design one needs a powerful workflow engine. This is crucial
was set to 0.1. The Powell method was used for a local search. The for combining different building blocks in Drug Design in a flexi-
number of initial points generated by Monte-Carlo methods was ble and configurable way. We plan to use P-PGRADE portal (Far-
set to 64. kas 2011) which is capable to construct complex workflows and
harness different types of grids and clouds.
Table 5. Testing results for Monte-Carlo method coupled with MBH me-
thod 7 CONCLUSIONS
The docking problem is basically a global optimization problem
Name Avrg. Energy Best Energy Avrg. Time (s) (1) where objective f x is a protein-ligand interaction energy.
Possessing the powerful techniques for resolving this problem is
Monte-Carlo -138,66 -154,47 10913,7 crucial for the efficiency of docking. In this paper we considered
several optimization methods for problem (1). Numerical experi-
ments showed that the best results can be achieved by combining
The obtained results demonstrate that even a very simple globali-
some globalization techniques (e.g. Monte-Carlo method) and the
zation strategy gives a noticeable improvement over plain semi-
monotonic basin hopping semi-local method coupled with the
local and local methods. But the running time is considerably
Powell local search algorithm.
higher.
We also outlined the future software environment which will pro-
vide a consistent and convenient access for a wide range of re-
6 SCIENCE GATEWAY INTEGRATION searchers to the drug-design experimetns.
The ultimate goal of our work is to create software environment
for molecular simulation. Such environments are useful to isolate ACKNOWLEDGEMENTS
the end-user from technical details of the application.
This study was supported in part by grant 10-07-00595 from the
Russian Foundation for Basic Research and the DEGISCO FP/7
The use of high-performance computing in docking is inevitable
European (contract no. 261561).
because in practice millions of ligands have to be processed inde-
pendently. Docking is perfectly suitable for the desktop grid com-
puting (Kiss, Greenwel 2010, Kiss 2010). We are going to develop REFERENCES
this application and deploy it at ISA RAS desktop grid Kitchen, D.B.; Decornez, H.; Furr, J.R.; Bajorath, J. (2004) Docking and scoring in
virtual screening for drug discovery: methods and applications. Nat. Rev. Drug
(boinc.isa.ru/dcsdg) and on a combined infrastructure that connects Discov., 3, 935-949.
ISA RAS desktop grid and service grid infrastructure provided by Zoete, V.; Grosdidier, A.; Michielin, O. (2009) Docking, virtual high throughput
EDGeS VO (www.edges.eu). Such combined infrastructures get screening and in silico fragment-based drug design. J. Cell. Mol. Med., 13, 238–
much attention in the grid community in recent years (Urbach 248.
Rosin, C.; Halliday, S.; Hart, W.; Belew, R. (1997) A comparison of global and local
2009). This approach provides the transparent seamless integration
search methods in drug docking. Proc 7th Intl Conf on Genetic Algorithms. San
of desktop and service grids and results in huge consolidated com- Francisco, CA., 221-228.
puting power. Morris, G.M.; Goodsell, D.S.; Halliday, R.S.; Huey, R.; Hart, W.E.; Belew, R.K.;
Olson, A.J. (1998) Automated docking using a Lamarckian genetic algorithm and
The deployment will done on top of the BNB-Grid (Evtushenko empirical binding free energy function. J. Comp. Chem., 19, 1639-1662.
2009) programming environment – a library for solving large-scale
3rd International Workshop on Science Gateways for Life Sciences (IWSG 2011), 8-10 JUNE 2011
Tavares, J.; Melab, N.; Talbi, E. (2008) An Empirical Study on the Influence of
Genetic Operators for Molecular Docking Optimization. Inria Research Report,
ISSN 0249-6399 ISRN INRIA/RR--6660--FR+ENG.
Tavares, J.; Melab, N.; Talbi, E. (2009) On the Efficiency of Local Search Methods
for the Molecular Docking Problem. EvoBIO. 104-115.
Solis, F.; Wets, R. (1981) Minimization by Random Search Techniques. Mathematical
Operations Research, 6: 19-30.
Nocedal, J.; Wright, S.J. (2006), Numerical Optimization (2nd ed.), Berlin, New
York: Springer-Verlag, ISBN 978-0-387-30303-1.
Leary, R. H. (2000) Global Optimization on Funneling Landscapes. J. of Global
Optimization 18, 4, 367-383.
Cole, J.C.; Nissink, J.W.M.; Taylor, R. (2005) Virtual Screening in Drug Discovery;
Alvarez, J.; Shoichet, B. Eds.; CRC Press: Taylor & Francis Group, Boca Raton,
London, New York, Singapore; pp. 379-415.
Morris, G.M.; Huey, R.; Lindstrom, W.; Sanner, M.F.; Belew, R.K.; Goodsell, D.S.;
Olson, A.J. (2009) AutoDock4 and AutoDockTools4: automated docking with se-
lective receptor flexibility. J. Comp. Chem., 30, 2785–2791.
Romanov, A.N.; Kondakova, O.A.; Grigoriev, F.V.; Sulimov, A.V.; Luschekina, S.V.;
Martynov Ya.B.; Sulimov V.B. ( 2008) The SOL docking package for computer-
aided drug design. Numerical Methods and Programming (Section 1. Numerical
methods and applications9(2), 64-84 (in Russian).
Romanov, A.N.; Jabin, S.N.; Martynov, Ya.B.; Sulimov, A.V.; Grigoriev, F.V.; Suli-
mov, V.B. ( 2004) Surface Generalized Born method: a simple, fast and precise
implicit solvent model beyond the Coulomb approximation. J. Phys. Chem. A,
108, 9323-9327.
Oferkin, I.V.; Sulimov, A.V.; Kondakova, O.A.; Sulimov V.B. (2011) Docking pro-
grams SOLGRID & SOL: parallel implementations. Numerical Methods and Pro-
gramming (Section 2. Programming) 12, 9-23 (in Russian).
Halgren, T.A. (1996) Merck Molecular Force Field. J. Comp. Chem., 17(5&6), 490-
641.
Ghosh, A.; Rapp, C.S.; Friesner, R.A. (1998) Generalized Born Model Based on a
Surface Integral Formulation. J Phys Chem B, 102, 10983-10990.
Polak, E. (1971) Computational Methods in Optimization. New York: Academic
Press.
Nash, S. (2000)A survey of truncated-Newton methods. Journal of Computational and
Applied Mathematics. Volume 124 Issue 1-2, Dec. 2000.
Powell M.J.D. (1964) An efficient method for finding the minimum of a function of
several variables without calculating derivatives, Computer Journal, 7 :155–162.
Wales, D.; Doye, J. (1997) Global Optimization by Basin-Hopping and the Lowest
Energy Structures of Lennard-Jones Clusters Containing up to 110 Atoms. J.
Phys. Chem. A, 101, 5111-5116.
Panteleev, A.V., Letova. T.A. (2005) Optimization methods via examples. Moscow,
544 p.
Yoo,M.S. et al. (2003) Oxidative stress regulated genes in nigral dopaminergic neur-
nol cell: correlation with the known pathology in Parkinson’s disease. Brain Res.
Mol. Brain Res., 110(Suppl. 1), 76–84.
Kim, J.B., Nhan, N.D. et al. (2008) DrugScreener-G: Towards an Integrated Environ-
ment for Grid-Enabled Large-Scale Virtual Screening and Drug Discovery. Pro-
ceedings ESCIENCE '08 Proceedings of the 2008 Fourth IEEE International Con-
ference on eScience, 666 – 671.
Kiss T., Greenwell P., Heindl H. et al. (2010) Parameter Sweep Workflows for Mod-
elling Carbohydrate Recognition, Journal of Grid Computing, volume 8, number
4, pp 587-601, DOI: 10.1007/s10723-010-9166-8
Kiss T., Szmetanko G. et al. (2010) Molecular Docking Simulations on a Combined
Desktop and Service Grid Infrastructure, Proceedings of the 3rd Grid Experience
workshop - Desktop Grid Applications for eScience and eBusiness, March 30,
2010, Almere, The Netherlands, p. 23-28
E. Urbach, P. Kacsuk, Z. Farkas et al. (2009) EDGeS: Bridging EGEE to BOINC and
Xtrem Web, Journal of Grid Computing, Vol 7, No. 3, pages 335 -354.
Evtushenko, Y., Posypkin, M., Sigal, I. (2009) A framework for parallel large-scale
global optimization. Computer Science - Research and Development 23(3), 211-
215.
Farkas, Z., Kacsuk, P. (2011) P-GRADE Portal: A generic workflow system to sup-
port user communities. Future Gener. Comput. Syst. 27, 5, 454-465.
DOI=10.1016/j.future.2010.12.001.
Evtushenko Yu. (1971) Numerical methods for finding global extreme (case of a non-
uniform mesh). U.S.S.R. Comput. Maths. Math. Phys., Vol. 11, pp. 1390-1403.
De Jong K. A. (2006) Evolutionary computation: a unified approach. MIT Press,
Cambridge MA