Merlitz, H.

Wenzel, W.

English

KITopen

Comparison of stochastic optimization methods forreceptor ligand dockingH. Merlitz, W. Wenzel *Forschungszentrum Karlsruhe GmbH, Institut fur Nanotechnologie, Postfach 3640, D 76021 Karlsruhe, Germany AbstractWe compare the efficiency of three stochastic optimization methods, simulated annealing, parallel tempering andstochastic tunneling to locate the global minima of complex and rugged potential energy surfaces arising from atomisticmodels for receptor ligand docking. The stochastic tunneling method proves to be the most efficient generic approachfor atomistic receptor ligand docking in the rigid ligand rigid receptor approximation.  2002 Elsevier Science B.V.All rights reserved.1. IntroductionThe development of methods to efficiently de-termine the global minima of complex and ruggedenergy landscapes remains a challenging problemwith applications in many scientific and techno-logical areas.In particular for NP-hard [1,2] problems,straightforward enumerative as well as sophisti-cated branch-and-bound techniques become pro-hibitively expensive and stochastic methods offerthe only acceptable compromise between thecomputational cost of the method and its reli-ability. In such techniques the global minimizationis accomplished by the simulation of a fictitiousdynamical process for a ‘particle’ on the multi-di-mensional potential energy surface (PES). Thisparticle explores the potential energy surface in abiased random walk that is designed to guide it tolow-energy regions.In one of the most widely used methods, thesimulated annealing technique (SA) [3], the PES isexplored in a Monte Carlo simulation with grad-ually decreasing temperature. The final tempera-ture is chosen small enough to constrain thedynamics of the particle to the immediate vicinityof the nearest local minimum of the PES, while thelargest temperature must be sufficiently high toallow an essentially random search of the PES. Onrugged PES this strategy routinely fails, becausetransition states between adjacent local minima aretoo high to be overcome at temperatures suffi-ciently small to resolve the energy differences be-tween them.Here we investigate a family of methods thathave sought to address this generic deficiency ofSA by allowing the dynamical process to passthrough thermodynamically inaccessible regions of* Corresponding author. Fax: +72 47 82 6434.E mail address: wenzel@int.fzk.de (W. Wenzel).the PES [4]. Generalized ensemble techniques suchas simulated tempering [5,6], parallel tempering(PT) [7] and multi-canonical algorithms [8,9], havebeen proposed to address this problem. In a pre-vious report, we have introduced the stochastictunneling method [10] (STUN) as an alternateapproach that avoids both the numerical cost ofmultiple simulations as well as the estimation ofauxiliary control fields. STUN was successfullyapplied to peptide folding [11], as were PT [12] orSA [13]. Here we compare this method with PTand SA in application to the receptor liganddocking (RLD) problem in an atomistic model.We find STUN to be an efficient and reliabledocking procedure, that succeeded even in caseswhere neither PT nor SA were able to locate theglobal minimum in repeated simulations.2. Receptor–ligand dockingIn the RLD problem, suitable ligands must beselected for a given, structurally characterized re-ceptor [14,15]. One approach is to screen largechemical databases in-silico for suitable ligands.For each ligand the best possible fit between ligandand receptor must be determined. Even in the mostsimple atomistic model, where both protein andligand are treated as inflexible molecules, efficientnumerical techniques to screen large databases inany reasonable timeframe are still lacking. Thereason for this difficulty lies in the competitionbetween two vastly different energy scales in theproblem, where steric repulsion competes withattractive electrostatic forces and hydrogen bond-ing to determine the global minimum of the PES.The tight fit between receptor and ligand (keylock-principle) complicates the optimizationproblem significantly because it is almost impos-sible to reorient the ligand within the receptor,while there are few specific interactions betweenligand and receptor outside the receptor pocket.For the simulations discussed below we used ascoring function:S ¼XProteinXLigandRijr12ij Aijr6ijþ qiqjerij!; ð1Þwhich contains the empirical Pauli repulsion (firstterm), the Van de Waals attraction (second term)and the electrostatic Coulomb potential (thirdterm). Neither solvation effects nor dielectricscreening were used in the simulations becausesuch terms alter the specifics of the affinity of agiven ligand to the receptor, but not the nature ofthe optimization problem. We note, however, thatthe inclusion of such terms generally leads tonondifferentiable potentials, which are not ame-nable to optimization methods such as basinhopping [16]. The ligands are simulated as rigidbodies, there are six degrees of freedom in thesimulations. In cases where rotatable bonds exist,the X-ray crystallographic structures of the dockedligands were taken from the PDB database. Theforce field parameters Rij and Aij are taken fromthe OPLSAA force field [17] and the scoringfunction is pre-calculated on grids. The atomicaffinity grids are interpolated using a logarithmicinterpolation technique [18].3. MethodsIn SA the PES is explored in a stochastic pro-cess, where a given configuration with energy E1 ismodified and the new configuration with energy E2is accepted with probability:p ¼ expðbðE1  E2ÞÞ; ð2Þwhere b ¼ 1=kbT is an inverse fictious temperature.For fixed b the simulation samples the configura-tion space with a thermodynamic equilibriumdistribution corresponding to this temperature. InSA, the temperature is gradually reduced to zeroaccording to a cooling schedule, so that the systemshould end in principle in its ground state, orglobal energy minimum at the end of the simula-tion. For rugged or glassy PES, however, ergod-icity is routinely lost in the cooling process, which‘freezes’ in some metastable state.For such PES, the freezing problem may becircumvented by allowing a trapped particle toescape from a local minimum by increasing thetemperature of its simulation. Following this ideathe parallel tempering method replaces the unidi-rectional cooling of SA by a set of concurrentH. Merlitz, W. Wenzel4. ResultsRLD simulations were carried out for threerepresentative natural ligands, benzamidine, met-hotrexate and retinol in their natural receptors, b-trypsin (3 ptb), dihydrofolate reductase (4 dfr),and retinol binding protein (1 rbp), respectively. Ingeneral RLD simulations are complicated by thefact that (a) not much is known a priori about thecenter of mass location of the ligand with respectto the receptor pocket and (b) a large fraction ofthe search space is sterically forbidden.To succeed in this scenario we optimized thedynamics of the underlying random walks as fol-lows: the ligand is embedded into a symmetricalellipsoid, where analytical expressions for thefriction coefficients ft (translation), fa (rotationaround long axis a) and fb (rotation around shortaxis b) exist [20]. The Brownian dynamics tech-nique in the diffusive regime is based on the updatescheme [21]rðt þ dtÞ ¼ rðtÞ þ FðtÞdt=ft þ X; ð5Þwhere F is a drift term and X is the random dis-placement sampled from a Gaussian with zeromean and width hX 2i ¼ 2kBTdt=ft. Analogousequations exist for the rotational updates dxa anddxb around axes a and b, respectively. Takingretinol, the ratio a=b  4:5 is reflected in the dif-ferent friction coefficients fi which determine thevalues for the angular shifts, leading to a ratiodxa=dxb  2:4. Similarly, the translational updateis physically related to the rotational updates andadapted to the specific shape of the ligand. Incontrast to Brownian dynamics, the time step dthas no physical counterpart but can be used as afree parameter to adjust the acceptance rate.For the drift term we introduce a point psomewhere inside the cavity of the receptor. Thedrift force dfd  kdT ðr pÞ defines an addi-tional, systematic contribution to the dynamics.The strength of the drift is proportional to thedistance jr pj as it were in a harmonic oscillatorfield. If there were no further external forces, thisdrift would lead to a Gaussian localization of thecenter of mass coordinates of the ligand [22]. Theadvantage of this approach over the introductionof a penalty function is that the structure of thepotential surface remains unchanged. No morethan a bias in the sampling procedure is added tothe random displacements. The ligand is localizedin a way that its probability distribution fills thecavity of the receptor and is significantly reducedfar outside the cavity. The localization volumedoes not depend on the step size dt and the tem-perature T.In all simulations ligands were placed in arandom position outside the cavity and we aver-aged the results of 50 runs of predescribed stepnumber. A ligand was defined as ‘docked’ if theaverage RMS deviation of the atoms from theglobal minimum was less than 0.1 nm. The po-tential values were pre-calculated on cubic gridswith a grid constant of 0.04 nm and a dimension of3  3  3 nm3. Generally, a grid whose size is welladapted to the cavity dimension can improve theefficiency of the simulation, but in the presence ofthe drift term the actual size and shape of the gridbecomes secondary. As a drift center we havechosen the actual center of mass of the docked li-gand, but due to the diffuseness of the localizationvolume there is no need for an accurate estimatefor that. Any drift term which is centered near thecavity and whose delocalization allows the ligandto probe the entire cavity allows for the desireddocking.The optimum cooling schedule for SA was tostart with an initial temperature of 3000 K andthen cooling down to 300 K toward the end of therun. In PT, we face the problem that n simulta-neous simulations produce an n-fold increase ofthe number of steps for a given simulation length.A certain minimum simulation length, however, ismandatory for the ligands to find their way intothe cavity. This imposes the upper limit of threeconfigurations for PT to be competitive withSTUN and SA, a number too low to exploit thefull potential of this technique. The temperatureswere found to be best at 1000, 1500 and 2250 K forthe three configurations.4.1. Benzamidine dockingBenzamidine is an oblate shaped ligand which,in the ellipsoidal approximation, has an axial ratioof a=b ¼ 11, i.e., it is a flat discus. The docking toH. Merlitz, W. Wenzel

Comparison of stochastic optimization methods for receptor-ligand docking

https://publikationen.bibliothek.kit.edu/110052725/149885753

Comparison of stochastic optimization methods for receptor-ligand docking

Abstract

Similar works

Full text

Available Versions

KITopen