Solving molecular flexible docking problems with metaheuristics: A comparative study

Solving molecular flexible docking problems with metaheuristics: A comparative study Esteban López-Camacho, María Jesús García Godoy, José García-Nieto, Antonio J. Nebro , José F. Aldana-Montes Lenguajes y Ciencias de la Computación, Universidad de Málaga, Bulevar Louis Pasteur 35, Málaga, Spain Keywords: Molecular docking Optimization Metaheuristics Experimental comparison a b s t r a c t The main objective of the molecular docking problem is to find a conformation between a small molecule (ligand) and a receptor molecule with minimum binding energy. The quality of the docking score depends on two factors: the scoring function and the search method being used to find the lowest binding energy solution. In this context, AutoDock 4.2 is a popular C++ software package in the bioinformatics community providing both elements, including two genetic algorithms, one of them endowed with a local search strategy. This paper principally focuses on the search techniques for solving the docking problem. In using the AutoDock 4.2 scoring function, the approach in this study is twofold. On the one hand, a number of four metaheuristic techniques are analyzed within an extensive set of docking problems, looking for the best technique according to the quality of the binding energy solutions. These techniques are thoroughly evaluated and also compared with popular well-known docking algorithms in AutoDock 4.2. The metaheuristics selected are: generational and a steady-state Genetic Algorithm, Differential Evolution, and Particle Swarm Optimization. On the other hand, a C++ version of the jMetal optimization framework has been integrated inside AutoDock 4.2, so that all the algorithms included in jMetal are readily available to solve docking problems. The experiments reveal that Differential Evolution obtains the best overall results, even outperforming other existing algorithms specifically designed for molecular docking. 1. Introduction Molecular docking tools are used to predict the conformation and relative orientation of two or more molecular constituents. Pioneered in the early 1980s, these computational approaches have made an important contribution to the location and characterization of binding sites in macromolecules, as well as in the discovery, design and comparison of new drugs to find new compounds with high affinity to the targeted proteins [1]. The objective of docking is to obtain a model of interaction between a ligand and a macromolecule, characterized by a minimum binding energy. This model is based on known ligand and macromolecule 3D structures. In the past, this requirement was a problem due to the lack of 3D crystal structures. However, with the development of techniques such as high-throughput protein ∗Corresponding authors. Tel.: +34 952133310. purification, the X-ray crystallography and nuclear magnetic resonance (NMR) spectroscopy, the number of structures stored in the Protein Data Bank database (PDB) [2] has increased. This database currently contains 100,547 entries for 3D structures obtained by X-ray crystallography (over 80%), NMR (about 16 %) and theoretical modeling (2%). Most of these PDB structures have an important function in the metabolic and biosignaling pathways and can be considered as drug targets in docking analysis. In addition, the advances in human genome sequencing have provided a basis for the detection of 3D protein structures related to the regulation and the expression of genes. In this context, these 3D structures offer a high diversity of biologically active macromolecules which are used to generate drug interaction models for the discovery of new targets. In the computational development of docking software, researchers in this field have traditionally focused on two of the components which determine the quality of the docking software: the scoring function and the optimization algorithm. The scoring function enables the evaluation of conformations between ligand and protein, to detect those conformations with the minimum binding energy. Moreover, the conformational model used in E-mail addresses: [email protected] (E. López-Camacho), [email protected] (M.J. García Godoy), [email protected] (J. García-Nieto), [email protected] (A.J. Nebro), [email protected] (J.F. Aldana-Montes). docking software tools allows the rotation around torsional degrees of freedom of the flexible ligand and side chains of receptor in order to create more realistic docking simulations [3]. So, given a scoring function, the goal is to use an optimization algorithm to find the best docking solutions which are those complexes with the minimum binding energy. In this paper, we focus on metaheuristics [4], a broad family of non-exact optimization techniques, including Evolutionary Algorithms (EAs), Particle Swarm Optimization (PSO) [5], Ant Colony Optimization (ACO) [6], Simulated Annealing (SA), Tabu Search (TS), and many others. Genetic Algorithms (GAs), a subfamily of EAs, are the most well-known and used metaheuristics. A popular tool for molecular docking is AutoDock [7], a C++ software package which in its 4.2 release provides a SA and two GAs algorithms, one of which, referred to as the Lamarckian GA, incorporates a local search [8]. These techniques involve specific operators using problem information, which leads them to outstanding results for rigid and flexible molecular confirmations. The main motivation in this study is to analyze the performance of a set of general-purpose metaheuristics (in their canonical design) to determine whether any of them can lead to better scoring values compared to the techniques already provided by AutoDock. In particular, we focus on two variants of standard GA (generational and steady-state), a PSO algorithm, and a Differential Evolution (DE) solver [9]. To achieve this goal, instead of trying to incorporate the new algorithms into the source code of AutoDock, the approach has been, to use a software framework oriented to solving optimization problems with metaheuristics. Specifically, the framework used is jMetal [10], a Java-based object-oriented software library aimed at multi-objective optimization but it also incorporates a number of single-objective algorithms. As AutoDock and jMetal are written in different programming languages, to obtain an efficient implementation we have opted to develop a C++ version of jMetal. The result is a software package that integrates both tools [11], in such a way that the algorithms available in jMetal can be used to optimize the binding energy function of AutoDock 4.2. Consequently, AutoDock users also benefit because they have the choice of using different optimization techniques beyond the ones already incorporated inside the software, and researchers in metaheuristics can use a real world problem, like molecular docking in their work. For the algorithmic evaluation and comparison, a thorough experimentation is performed in this study, where algorithms are subjected to a benchmark of molecules with different properties of flexibility, size and resolution. As mentioned, the flexibility in macromolecule and ligands is considered a more realistic docking problem and also a challenge given the more complexity of the problem, therefore, we have chosen to focus here on solving this problem. The ligands are of different sizes from small to large including cyclic urea ligands. The main contributions of this work can be summarized as follows: •Four metaheuristics (two GAs, PSO, and DE) are used to solve moleculardockingproblemsandtheirbehavioranalyzedinterms of performance. In addition, these algorithms are also compared with two well-known metaheuristics included in AutoDock 4.2 (GA and LGA). •A set of 83 docking instances is generated based on real existing molecules, then constituting an extensive benchmark for assessing algorithms. From this set, 75 PDB structures (instances) are used for performance assessment involving flexibility in HIV- protease macromolecules and ligands, whereas a subset of 11 molecular instances are used for the parameter setting (three instances are used in both subsets), revealing new configurations for fine-tuning algorithms in the context of flexible molecular docking. The benchmark is online available and adapted to be reproduced using compatible software.1Along this paper, all necessary information is provided to make this study replicable according to standard guidelines [12]. •A thorough experimentation with statistical assessment of results is performed to compare algorithms. An additional analysis in terms of convergence behavior is also provided. •A series of relevant solutions are analyzed from the point of view of their biological meaning, so as to offer insightful results for experts in the domain of application of molecular modeling and design. The remainder of this paper is organized as follows. Section 2 includesareviewofthestateoftheartconcerningotherapproaches to the problem being tackled here. Section 3presents the molecular docking problem, with solution encoding and fitness energy function formulation. In Section 4, the metaheuristics evaluated in this study are described. Section 5includes the description of the strategy followed, based on the jMetal + AutoDock integration. Section 6 reports the experimental setup, with parameter settings and problem instances. Results, comparisons, and analysis are presented in Section 7. Section 8concludes the paper and presents an overview of future work. 2. Related work Over the last two decades, different metaheuristics have been applied as search methods to solve the docking problem [13].One example is the docking software AutoDock, which incorporates three metaheuristic techniques. AutoDock is considered to be the most cited and one of the most used software packages in molecular modeling studies [7], because it is an efficient tool for virtual drug screening [14]. AutoDock was released in 1990, and it included a rapid search method using Monte Carlo simulated annealing [15]. However, this method proved to be inadequate for ligands with more than eight rotatable bonds [8]. Eight years later, in an attempt to improve the software, AutoDock 3.0 was released, adding the Genetic Algorithm (GA) and the Lamarckian Genetic Algorithm (LGA), which incorporates a local search, and an empirical binding free energy force that enables the prediction of binding of free energies. Docking analyses have demonstrated that the LGA is the most efficient search method of the three AutoDock algorithms in terms of the lowest energy found in a number of energy function evaluations [8]. AutoDock 4 was presented in 2009 [3]. It allows conformational models of side chains of proteins, provides torsional degrees of freedom and tries to solve the problem of flexibility in the receptor, a challenge in docking approaches. More recently, a new release has appeared, AutoDock 4.2, which incorporates several enhancements over AutoDock 4. The latest version includes a default unbonded state, different to the extended unbonded state of AutoDock 4, an improvement over the time required to run a high-quality docking with flexible and rigid components, involving an attempt to ensure compatibility between the different releases of AutoDock software. In 2010, as an improvement on the previous releases, the AutoDock authors implemented a new program for molecular docking called AutoDock Vina [16], containing differences with respect to the previous versions, such as the energy function, multithreading to speed up the execution in multicore processors, and the use of an iterated local search algorithm (ILS) as the search engine. The stopping condition of the ILS is adaptively determined, thus making it difficult to compare it fairly with other techniques that use a fixed number of function evaluations. 1http://khaos.uma.es/AutodockjMetal/instances.jsp. A number of approaches can be found in the current literature that have proposed metaheuristic techniques designed around AutoDock versions. Atilgan et al. [17] developed a new program named AutoDockX which incorporates a sustainable GA, namely Age-Layered Population Structure (ALPS), including the age attribute for individuals. Chen et al. [18] presented an algorithm called SODOCK, which is an adaptation of PSO including Solis and Wets local search and uses the energy function of AutoDock 3.05. Two other PSO related proposals are the varCPSO-ls algorithm, an extension of the CPSO algorithm with a local optimizer which is embedded inside the AutoDock 3 source code and uses its energy function [19], and the FIPSDock algorithm, which adopts the AutoDock 4.2 energy function [20]. DE has also been applied in this context. A first attempt is DockDE [21], a variant of DE which uses a older version of the AutoDock energy function. ODE is also an extension of DE enhanced by a local search algorithm and a pseudoelitism operator, and using the AutoDock scoring function [22]. A more recent version is SADock [23], that incorporates a Hooke Jeeves local search. Among studies on molecular docking with metaheuristics not based on AutoDock, there are several approaches that are also worth mentioning. PARADockS is a framework implemented to predict the ligand–protein interaction adapting PSO to several objective functions [24]. An Ant Colony Optimization approach is also presented in [25] using a systematic molecular simulator. A variant of DE called MolDock was parallelized on both GPU and CPU using a fitness function designed by the authors [26]. However, although the technique is adapted to a flexible docking receptor, it has not been evaluated using flexible targets. Other optimization methods such as multi-scale optimization models and information entropy-based searching techniques with narrowing space were applied in a new docking algorithm [27]. Herberlé et al. [28] review EAs applied to Mycobacterium tuberculosis docking targets. In terms of analysing the influence of algorithm operators and parameters, a study developed by Thomsen [21] compared the performance of the LGA and DockEA algorithms by selecting different EA operators, populations and usage of local search. However, the parameter setting study proposed, included very few docking problems and was only applied to the DockEA algorithm. Another interesting study performed by Tavares et al. [29] investigated the effects of Gaussian and Cauchy mutation operators through a locality analysis (small genotype variations imply small variations in phenotype); the results showed that Gaussian-based operators had a stronger locality than Cauchy-based operators. They also demonstrated that the results of runs using the Gaussian-based operator were better than those returned by the Cauchy-based operator. Some studies have been developed over the last five years based on multi-objective optimization, such as the reviews on the multiobjective application in different fields proposed by Nicolaou et al. [30] and Nicolotti et al. [31]; a new hybrid algorithm proposed by Grosdidier et al. [32] which uses two fitness functions where two variables are optimized; Sandoval-Perez et al. [33] used the NSGA-II algorithm implemented in the jMetal framework in order to optimize two variables (bonded and non-bonded energy terms); Oduguwa et al. [34] applied the NSGA-II, PAES, SPEA evolutionary multi-objective techniques optimizing up to three objectives of the energy function; and finally, Janson et al. [35] applied a multiobjective optimization approach using the PSO algorithm. All these approaches have addressed different aspects of the molecular docking optimization. However, three common features can be found in most of them which clearly differentiates them from the main focus in this paper: •Our study focuses on carrying out a comparative analysis with canonical metaheuristics in the solution of molecular docking problem. To the best of our knowledge, there have been no other studies performed which carry out this comparative analysis. •Most of the work presented in this section uses a set with a low number of structures to be evaluated and concentrates on rigid macromolecules to be docked. The current study includes an extensive set of 83 molecular structures. In this set, receptors are HIV-proteases and ligands that have different sizes and flexibility applied to macromolecule and ligand. These problems of macromolecule flexibility are not considered in the current reviewed literature, therefore we are providing new fresh results for the scientific community. •A systematic parameter setting with multiple combinations of parameter values statistically validated is performed here. In addition to all these contributions, a software framework has been developed,2which integrates AutoDock 4.2 and jMetal, enabling the community to use those algorithms that are not available in AutoDock. 3. The problem: molecular docking The main objective in a molecular docking problem is to find an optimized conformation between the ligand (L), a rigid/flexible smallmolecule, and the receptor (R), a rigid/flexible macromolecule that results in a minimum binding energy. The interaction between ligand and receptor can be described by an objective function calculated according to three components representing degrees of freedom: (1) the translation of the ligand molecule, involving the threeaxisvalues(x,y,z)inCartesiancoordinatespace;(2)theligand orientation, modeled as a four variables quaternion including the angle slope (w); and (3) the flexibilities, represented by the free rotation of torsion (dihedral angles) of the ligand and sidechains of the receptor. Solution encoding. As illustrated in Fig. 1, each problem solution for AutoDock and jMetal is encoded by a real-value vector of n+ 7 variables, in which the first three values correspond to the ligand translation, the next four values correspond to the ligand and/or macromolecule orientation, and the remaining nvalues are the ligand torsion dihedral angles. Nevertheless, with the aim of reducing the computational cost, a grid-based methodology has been implemented where the protein active site is embedded in a 3D rectangular grid, and on each point of the grid, the electrostatic interaction energy and the van der Waals terms for each ligand atom type are pre-computed and stored, taking into account all the protein atoms. In this way, the protein contribution at any given point is obtained by tri-linear interpolation in each grid cell. This interpolation leads to a range of translation variables (x,y,z) of 120 grid spacing points dimension (see [36]). The values of these variables are delimited between the range of the coordinates of the grid space that has been chosen for each problem. All ranges are selected randomly, so if the center of the grid is for example the (10, 10, 10) point, a solution with values of ten in its x,yand zwill be in such a position. In the case of orientation (quaternion) and torsion variables, they are measured in radians and encoded in the range of [−,]. Fitness function. In this approach, the energy scoring function (or fitness function) recommended in AutoDock 4.2 [36] is used to measure the quality of binding in solutions, as it was empirically determined in the scope of this tool. This scoring function, also called semi-empirical energy force field, computes the total free binding energy of the ligand-macromolecule complex according to the following equations: 2Available at URL: http://khaos.uma.es/AutodockjMetal/. Fig. 1. Solution encoding in AutoDock 4.2 and jMetal. The first three values (translation) are the coordinates of the center of rotation of the ligand. The next four values (quaternion) are the unit vector describing the direction of rigid body rotation (x,yand z) and the rotation of the angle degrees (w) that are applied. The rest of the values hold the torsion angles in degrees, being nthe number of torsions of the ligand. G =(VL−L bonded −VL−L unbonded)+(VR−R bonded −VR−R unbonded) +(VR−L bonded −VR−L unbonded +Gconf) (1) V=Wvdw i,j (Aij r12 ij −Bij r6 ij )+Whbond i,j E(t)Cij r12 ij −Dij r10 ij  + + Welec i,j qiqj ε(rij)rij +Wsol i,j (SiVj+SjVi)e(−r2 ij /22)(2) Gconf =Wconf Ntors (3) As formulated in Eq. (1), the free binding energy function is calculated from the differences between ligand/s (L) and receptor/s (R) in bonded and unbonded states. That is, the free energy of binding is based on the evaluation of the transition of the ligand and protein intramolecular energetics from an unbonded to a bonded state, and the intermolecular energetics of ligand–protein complex. Therefore, the force field involves six pair-wise evaluations (V) plus a term of conformational entropy (Ntors), which is directly proportional to the number of rotatable bonds of the ligand molecule (Eq. (3)). Each pair of energetic evaluation terms includes evaluations of dispersion/repulsion (vdw), hydrogen bonds (hbond), electrostatics (elec), and desolvation (sol). Weights Wvdw,Whbond,Wconf, Welec, and Wsol of Eqs. (2) and (3), are constants for van der Waals, hydrogen bonds, torsional forces, electrostatic interactions, and desolvation, respectively. In Eq. (2),rij represents the interatomic distance, Aij and Bij in the first term are Lennard–Jones parameters taken from Amber force field [37]. Similarly, Cij and Dij in the second term are Lennard-Jones parameters for maximum well depth of potential energies between two atoms, and E(t) represents the angle-dependent directionality. The third term in Eq. (2) uses a Coulomb approach for electrostatics. Finally, the fourth term is calculated from the volume (V) of the atoms that are surrounding a given atom weighted by S, and an exponential term which involves atom distances. Ligand–protein docking is a highly complex optimization problem, with unknown optimum and usually characterized by multimodal landscape energy functions [38]. In addition, the computational cost of each energy evaluation increases with the number of atoms in complex ligand-protein (with thousands of them), hence involving millions of energy evaluations, since a minimum quality of molecular binding is mandatory in chemistry research. Therefore, the use of metaheuristic approaches is highly recommendable for molecular docking, since they are able to explore a great number of combinations with a fast convergence to successful solutions [4]. 4. Metaheuristics This section describes the four metaheuristic algorithms evaluated in this study. Specifically, they are two evolutionary algorithms, steady state Genetic Algorithm (ssGA) and generational Genetic Algorithm (gGA); and two differential vector techniques, Particle Swarm Optimization (PSO), and Differential Evolution (DE). These techniques have been selected as they are widely used in the literature to solve real-coded optimization problems [39]. In addition, they lead to experiments with different population structures and reproduction mechanisms. 4.1. Genetic Algorithms Genetic Algorithms [40] are the most popular metaheuristic algorithms belonging to the subfamily of Evolutionary Algorithms (EAs). A GA iterates a process in which two solutions (parents) are selected from the whole population with a given selection criterion. These parents are then recombined to obtain offspring solutions. The offsprings are mutated and finally they are evaluated and inserted back into the population following a given replacement criterion. As the recombination operator, a simulated binary crossover (SBX) is used for continuous variables [39]. The mutation process is carried out by selecting one of the elements in the solution, and assigning a new value in the range of problem variables. Algorithm 1 summarizes the operations of a canonical GA. Algorithm 1. Pseudocode of GA 1: P0←initializePopulation() 2: while g< MAXIMUMgdo 3: Rg←recombine(Pg) 4: Mg←mutate(Rg) 5: evaluate(Mg) 6: Pg+1 ←select(Mg∪Rg) 7: end while There are two main versions of GA: steady state GA (ssGA) and generational GA (gGA). The difference between the ssGA and the gGA is the way in which the population is being updated with the new individuals generated during the evolution. In the case of ssGA, the new individuals are directly inserted into the current population, whereas in the case of the gGA, a new auxiliary population is built with the offsprings obtained and then, once this auxiliary population is full, it replaces the current population. Thus, in ssGAs the population is being asynchronously updated with the newly generated individuals, while in the case of genGAs all the new individuals are updated at the same time, in a synchronous way. It is worth mentioning that jMetal [41] incorporates a number of GAs, including steady-state (ssGA), generational (gGA), as well as several cellular GA variants. In this approach, we have selected the ssGA and the gGA, which have been configured with common settings used in many approaches [39,42]: binary tournament as the selection operator, and polynomial mutation and simulated binary crossover (SBX). All these operators work without any specific knowledge of the problem. In addition, AutoDock 4.2 also incorporates a GA following a steady-state selection scheme, where the degree of elitism (i.e., the number of best solutions that survive into the next generation) is user-defined. The evolutionary operators are proportional selection, 2-point crossover, and Cauchy mutation. The crossover operators are adapted to the problem solution structure, as the solution is never split using a quaternion value position. Furthermore, a hybrid version of GA with local search, the aforementioned Lamarckian GA (LGA), is also provided in AutoDock 4.2. The local search method follows a pseudo-Solis–Wets scheme [43], and consists of a number of iterations during which rotational, translational, and torsional degrees of freedom of the ligand are randomly sampled at a given step. These last two GAs are also compared with the four selected metaheuristics in the experimental sections of this study. 4.2. Particle Swarm Optimization Particle Swarm Optimization [5] is a population based metaheuristic inspired by the social behavior of birds within a flock, and was initially designed for continuous optimization problems. In this algorithm, each particle position xi(solution) is updated each generation tby means of Eq. (4). xt+1 i=xt i+vt+1 i(4) where vt+1 iis the velocity vector of the particle given by vt+1 i=(vt i+Ut[0, ϕ1]·(pt i−xt i)+Ut[0, ϕ2]·(bt i−xt i)) (5) In this formula, pt iis the personal best position the particle ihas ever stored, bt iis the position found by the member of its neighborhood that has had the best performance so far. Acceleration coefficients ϕ1and ϕ2control the relative effect of the personal and social best particles, and Utis a diagonal matrix with elements distributed in the interval [0, ϕi], uniformly at random. Finally, is the Clerc’s constriction coefficient [44]. This coefficient is calculated by means of Eq. (6), from the two acceleration coefficients ϕ1and ϕ2, being the sum of these two coefficients that determines the to use. Usually, ϕ1=ϕ2= 2.05, giving as results ϕ= 4.1, and = 0.7298 [44] =2 |2−ϕ−ϕ2−4ϕ| , with ϕ = i ϕi, and ϕ > 4 (6) Algorithm 2. Pseudocode of PSO 1: S←initializeSwarm() 2: while t< MAXIMUMtdo 3: for each particle xt ido 4: updateVelocity(vt i) //Eq. (5) 5: updatePosition(xt i) // Eq. (4) 6: evaluate(xt i) 7: update(pt i) 8: end for 9: updateLeader(bt i) 10: end while Algorithm 2 describes the pseudo-code of PSO. In concrete, this PSO specification is based on the canonical version proposed by Cler’s [44], for which only one specific parameter is required to be tuned (ϕ). The algorithm starts by initializing the Swarm (Population) which includes both the positions and velocities of the particles. The corresponding piof each particle is randomly (uniformly) initialized, as well as the leader b. Then, during a maximum number of iterations, each particle flies through the search space updating its velocity and position (Lines 4 and 5), it is then evaluated (Line 6), and its piis also calculated (Line 7). At the end of each iteration, the leader bis updated. 4.3. Differential Evolution Differential Evolution [9] is a stochastic population based algorithm designed to solve optimization problems in continuous domains. The task of generating new individuals is performed by differential operators such as the differential mutation and crossover. A mutant individual wt+1 iis generated by the following Eq. (7): wt+1 i←xt r1+·(xt r2−xt r3)(7) where r1, r2, r3, r4∈{1, 2, . . .,i−1, i+ 1, . . .,N}are random integers mutually different, and also different from the index i. The mutation constant > 0 stands for the amplification of the difference between the individuals xrsand it avoids the stagnation of the search process. Each mutated individual undergoes a crossover operation with the target individual xt i, by means of which a trial individual ut+1 iis generated. A randomly chosen position is taken from the mutant individual to prevent that the trial individual replicating the target individual. ut+1 i(j)←wt+1 i(j) if r(j)≤Cr or j=jr, xt i(j) otherwise. (8) As shown in Eq. (8), the crossover operator randomly chooses a uniformly distributed integer value jrand a random real number r∈[0, 1], also uniformly distributed for each component jof the trial individual ut+1 i. Then, the crossover probability Cr, and rare compared just like jand jr. If ris lower or equal than Cr (or jis equal to jr) then we select the jth element of the mutant individual to be allocated in the jth element of the trial individual ut+1 i. Otherwise, the jth element of the target individual xt ibecomes the jth element of the trial individual. Finally, a selection operator decides on the acceptance of the trial individual for the next generation if and only if it yields a reduction (assuming minimization) in the value of the fitness function f(), as shown by the following Eq. (9): xt+1 i←ut+1 iif f(ut+1i) ≤f(xt i), xt iotherwise. (9) Algorithm 3 shows the pseudocode of DE. After initializing the population, the individuals evolve during a number of iterations (MAXIMUMt). Each individual is then mutated (Line 5) and recombined (Line 6). The new individual is selected (or not) following the operation of Eq. (9) (Lines 7 and 8). Algorithm 3. Pseudocode of DE 1: P←initializePopulation(Ps) 2: while t<MAXIMUMtdo 3: for each individual iof Pdo 4: choose mutually different rsvalues 5: wt+1 i←mutation(xt rs, ) //Eq. (7) 6: ut+1 i←crossover(xt i,wt+1 i, Cr) //Eq. (8) 7: evaluate(ut+1 i) 8: xt+1 i←selection(xt i,ut+1 i) //Equation 9 9: end for 10: end while 5. Optimization strategy: jMetal + AutoDock The proposed optimization strategy is composed of two main processes acting in parallel: an optimization algorithm and a molecular docking procedure. The optimization part is carried out through a metaheuristic algorithm (ssGA, gGA, PSO, or DE) provided by jMetal [41]. jMetal is an open-source framework intended for single/multi-objective optimization metaheuristics. The molecular binding procedure is performed by AutoDock 4.2 [36], a C++ tool for virtual drug discovery widely used for real cases which involves rigid, as well as flexible docking. The proposed strategy implements an integration between jMetal and AutoDock in such a way that the algorithms available in jMetal can communicate with AutoDock to cooperate on the molecular docking optimization. To carry out this integration, jMetal framework was translated to C++ from scratch. This new version in C++ was called jMetalCpp. Such a solution consists of running AutoDock and jMetalCpp in two different threads inside the same process, memory is therefore shared so they can communicate with each other and they synchronize using mutexes. This approach is efficient and flexible, allowing any of the algorithms included in jMetalCpp to be easily used for solving docking problems. In this way, as illustrated in Fig. 2, when a given metaheuristic algorithm in jMetal is executed, it performs the optimization procedure and generates new solutions. Whenever the binding quality of a new solution has to be numerically quantified, it is sent to AutoDock to be evaluated with the binding energy function (Eq. (1)). After this evaluation, AutoDock returns the corresponding binding energy to the algorithm so it can assign this fitness value to the evaluated solution. This process is repeated until reaching a preestablished stop condition. As a result, the best solution found so far by the metaheuristic is returned to AutoDock to generate the output data concerning the optimized docking. Thus, the final results follow the standard format tailored to AutoDock users. 6. Experimental setup This section describes the sets of molecular (ligand–receptor) problem instances used for optimal docking in this study. Then, the methodology for preparing docking instances and software packages is also explained. Finally, the parameter settings of the evaluated algorithms are described. 6.1. Molecular docking problem instances The experimental procedure followed in this study comprises two different sets of protein–ligand complexes. Both sets of molecules were taken from the PDB database [45], and they are available online for experimental reproduction. The first set consists of 11 molecules with flexible ligands and receptors. This set has been used to carry out the parameter setting of algorithms (DE, PSO, ssGA, and gGA) evaluated in multiple systematic docking simulations. The second set comprises 75 molecules, also with flexible ligands and macromolecules, and is used for experimental comparisons. Specifically, the first set of molecules is used to train the learning procedure induced by tuning parameters in algorithms, whereas the second set is used to test and validate the actual performance of selected metaheuristics on a great number of docking instances. Table 1 summarizes the first set of selected problems showing the X-ray crystal structures, the PDB accession code, the structure resolution, and the reference where they were initially analyzed. In concrete, this first set of instances (1AJV, 1AJX, 1D4K, 1G2K, 1HIV, 1HPX, 1HTF, 1HTG, 1HVH, 1VB9, and 2UPJ) has been used for parameter settings, since they constitute a varied number of complexes with heterogeneous features and sizes. In addition, these instances have been previously used in [3] where the macromolecule flexibility is incorporated in the energy function of AutoDock 4.2. According to the size and type of ligand, 1HPX, 1HIV, 1HTG and 1D4K complexes include large inhibitors (the largest being ligand 1HPX and the smallest, 1D4K); the ligand of 1HTF structure is categorized as a small ligand and 1AJX, 1AJV, 2UPJ, 1HVH and 1G2K complexes include cyclic urea inhibitors. In general, the ligands’ torsional degrees of freedom are limited to 10, selecting those torsions that allow the fewest number of atoms to move around the ligand core. For the flexibility in the receptor, previous docking studies [3] with this set have demonstrated that residue ARG-8 showed the highest number of bad atom contacts with the ligand. Therefore, the flexibility docking was carried out adding three degrees of freedom in ARG-8, in both A and B receptor chains. The second set consists of an extensive number of complexes (75) used by [3]. Similarly to the first set, the flexibility of ligand of this second set was limited to 10 torsional degrees of freedom, and the macromolecule flexibility is performed adding three degrees of freedom in ARG-8 residue in both side chains, A and B. This aminoacid was selected due to previous docking studies in which these complexes showed the highest number of bad atom contacts with the ligand with resolution ( ˚ A) ranges between 1.09 and 2.8 (more details in [3]). As shown in Table 2, this second group of complexes can be classified in terms of type of ligand, from small to large size inhibitors, and cyclic urea inhibitors. The set of structures classified according to the ligand size, the PDB code, and the range of crystalographic resolution in ˚ A are also shown in Table 2. In general, a total number of 83 structures have been used in experiments, since three instances in the set of 11 molecules are also used in the set of 75. Therefore, the parameter setting experiments steer the algorithms being compared to sources of knowledge, both dependent and independent. These sets of problems are complex and widely demanded by the docking user community, although it constitutes an extensive and varied enough set to test the compared techniques from an algorithmic point of view. 6.2. Carrying out the docking experiments Solving the docking experiments with AutoDock and jMetal requires a number of steps that are enumerated here: Fig. 2. This schema shows the integration between AutoDock and jMetalCpp. Solutions generated by jMetalCpp framework are evaluated by AutoDock and sent back to jMetalCpp. Once the stop condition has been met, the best solution is returned to the AutoDock code which in turn generates an output file with the results. These results can be visualized with any tool that is suitable for visualizing AutoDock format (DLG) like AutoDockTool (ADT). Table 1 X-ray crystal structure coordinates taken from the PDB database and used for parameter settings, for flexible receptors. Their accession codes from the PDB database, resolution and references where they were initially analyzed are also presented. Protein–ligand complexes PDB Resolution ( ˚ A) References HIV-1 protease/AHA006 1AJV 2 [46] HIV-1 protease/AHA001 1AJX 2 [46] HIV-1 protease/Macrocyclic peptidomimetic inhibitor 8 1D4K 1.85 [47] HIV-1 protease/AHA047 1G2K 1.95 [48] HIV-1 protease/U75875 1HIV 2 [49] HIV-1 protease/KNI-272 1HPX 2 [50] HIV-1 protease/GR126045 1HTF 2.20 [51] HIV-1 protease/GR137615 1HTG 2 [51] HIV-1 protease/Q8261 1HVH 1.80 [52] HIV-1 protease/˛-D-glucose 1VB9 2.20 [53] HIV-1 protease/U100313 2UPJ 3 [54] 1. Preparing the macromolecule and the ligand: The X-ray crystallographicstructuresweredownloadedfromthePDBdatabase[45]. Chimera UCSF software [55] was used to remove small molecules such as solvent molecules, non-interacting ions, waters, etc. from crystal structure. After this, the AutoDockTools (ADT) graphical interface suite [36] is used to prepare the macromolecule/ligand coordinate files. For the preparation of ligand PDBQT files, partial atomic charges and the AutoDock atom types (for each ligand atom) are computed and assigned. In the case of macromolecule PDBQT files, atom partial charges and hydrogens are added by using, respectively, Gasteiger and Babel methods [36]. The flexibility of the ligand is limited to 10 torsional degrees of freedom and six for the macromolecule. The rigid and flexible parts of HIV-proteases are specified and saved as corresponding rigid and flexible PDBQT files. 2. Running Autogrid: AutoDock simulations require a previous calculation of grid maps in order to reduce the acting area for ligand–receptor movements. This grid calculation means that optimization algorithms deal with a reduced search space, thus guiding the search to more promising areas. These maps are calculated through Autogrid [36], once the coordinates have been established (120 ˚ A x 120 ˚ A x 120 ˚ A and 0.375 ˚ A of grid spacing) and centered on the cocrystallized ligand binding site. Once the Autogrid files have been obtained, a docking parameter file is configured to run AutoDock. Table 2 X-ray crystal structure coordinates taken from the PDB database and used in docking experiments. They consist of 75 molecules with accession codes from the PDB database. The range of resolution ( ˚ A) of each subgroup is also shown in the last column. Ligand type PDB code Resolution ( ˚ A) Small size 1B6M, 1HTE, 1HPV, 1B6L, 1BDL, 1HEG, 1HIH, 1A9M, 1GNM, 3AID, 1HBV, 1GNN, 1HSG, 1AAQ, 1TCX, 1GNO, 1BDQ, 1BDR, 1SBG, 1KZK, 1Z1H, 1Z1R, 1Z1H 1.09–2.8 Medium size 5HVP, 1K6P, 1MUI, 1JLD, 1B6J, 1K6C, 1IZI, 4PHV, 1HEF, 1K6T, 1IZH, 1HXW, 1IIQ, 1HPS, 1B6P, 1B6K, 4HVP, 2BPX, 2BPV, 1D4L, 1K6V, 1MTR, 1D4K 1.75–2.8 Large size 1VIK, 1HIV, 1HVS, 8HVP, 1ODY, 1A94, 3TLH, 1HVI, 1HOS, 1VIJ, 9HVP, 7HVP, 1HVJ, 1HVK, 1HVL, 1HWR, 1HTG 1.8–2.8 Cyclic urea 7UPJ, 1QBT, 1BV7, 1QBR, 1HPO, 1MEU, 1G35, 1PRO, 1BWA, 1BWB, 1QBU, 1DMP, 1BV9, 1MES 1.8–2.5 3. Running AutoDock + jMetal: Using the configuration files created in the previous steps, AutoDock is executed for each docking instance (as explained in Section 5). Original metaheuristics provided by AutoDock (GA and LGA) are executed here by means of AutoDock 4.2. Evaluated algorithms included in the jMetal framework (gGA, ssGA, DE and PSO) are executed in the proposed optimization strategy jMetal + AutoDock. All algorithms used, including the four from jMetal and the two from AutoDock, have been executed in computers equipped with modern dual core processors, 2 GB RAM, and Windows 7 O.S. They operate under a Condor [56] middleware platform, managing a maximum number of 400 processors, that acts as a distributed task scheduler (each task dealing with one independent run). Finally, given that metaheuristics are stochastic techniques, 30 independent runs have been carried out per algorithm configuration and instance, in order to obtain statistical confidence in the context of the whole experimentation. 6.3. Parameter settings In this study, parameters specific to each algorithm have been systematically tuned, whereas those similar ones involved in the experimental procedure have been commonly set, in order to establish comparisons as fairly as possible. Concerning specific parameters, Table 3 shows the results obtained from executions with all parameter combinations considered for each algorithm, when solving the set of 11 molecular instances (Table 1) used for this purpose (training).3These results are then compared by means of a Friedman statistical test [57] the ranking values of which are shown in this table for each pair of parameters considered, and for each algorithm evaluated. In this way, Genetic Algorithms in jMetal (gGA and ssGA) apply binary tournament selection, simulated binary crossover (SBX) with crossover probability pc= 0.9, and polynomial mutation with mutation probability pm= 1/L(Lis the number of problem variables). These two parameters (pcand pm) have been set with their standard values as proposed in the literature [39,42], leading to the systematic tuning for distribution indexes of crossover and mutation, cand m, respectively. Therefore, according to Table 3 (top left and right), distribution indexes resulted in values c= 5 and m= 100, in the case of gGA, and c= 5 and m= 5, in the case of ssGA. An interesting observation in this regard is that, low distribution indexes in SBX show the best performance for the two GAs in jMetal, therefore leading to an exploding behavior. In contrast, for gGA, the mutation distribution index is higher (m= 100), which induces a high diversity in the population, and then, an explorative search procedure. In the case of DE (variant rand/1/bin), mutation constant and crossover probability Cr are tuned as shown in Table 3 (bottom left), with the result that the combination of = 0.5 and Cr = 0.9 are the best performing parameters for this algorithm (Friedman’s rank 2.59). Finally, for PSO, acceleration coefficients (ϕ1and ϕ2) are tuned using all combinations of considered values (2.05, 2.10, 2.15, 2.20, and 2.25), performing a series of 30 independent runs for each one of the combinations. Therefore, as shown in Table 3 (bottom right), the best parameter combination found for PSO is ϕ1= 2.05 and ϕ2= 2.05, resulting in the lowest Friedman value (2.40) in this ranking. In addition, these results agree with the ones analytically set in the study of Clerc [44] concerning the PSO constriction factor. Another interesting observation in Table 3 lies in the resulting p-values of statistical tests for tuned algorithm, which are lower 3In total, these tuning experiments consisted of: 4 (algorithms) ×25 (parameter combinations) ×30 (independent runs) ×11 (problem instances) = 33,000 tuning experiments. than the confidence level (˛= 0.05) in almost all cases (excepting ssGA, although showing a low p-value = 2.01E−1). Therefore, it is notable that selected combinations of values are diverse enough to cover a wide range of possibilities in the parameter setting, hence obtaining fine-tuned parameters. Algorithms GA and LGA included in AutoDock 4.2 were set as done in the reference study in which they were proposed [8]. Both GAs apply 2-point crossover with a crossover probability pc= 0.8 and Cauchy mutation with probability pm= 0.02, ˛= 0.0, and ˇ= 1.0. The local search used in the LGA follows a pseudo Solis–Wets scheme with a maximum number of iterations of 300. Some parameters that are applied to both the algorithms in jMetal and AutoDock 4.2, such as the probabilities of crossover and mutation in the GAs, differin their values (e.g., pc= 0.9 injMetal andpc= 0.8inAutoDock), in order to maintain the default values adopted in each framework instead of unifying them. Table 4 summarizes the set of parameters selected to evaluate the algorithms. In this table, the population size is also included, which has been set to 150 individuals for all algorithms and instances. This population size is used as the default value in AutoDock [36], as in previous studies [58] where the same HIV- protease structures were tested. The number of energy function evaluations is 1,500,000 for all algorithms (including those from AutoDock4 docking simulations), then establishing a common stop condition for the experiments. In this regard, in the case of the instances of docking problems with flexible residues it is recommended that they be simulated using a large enough number of evaluations [3], so this value has been selected from preliminary trace observations (as shown in Fig. 4 and explained in next section), to assure convergence on a candidate solution, but avoiding taking an excessive time to do it. 7. Results In this section, the performance behavior of the metaheuristics studied is analyzed, including numerical and statistical comparisons. In addition, a biological analysis of selected outstanding solutions is also provided for the experts in the problem domain. 7.1. Performance comparisons Tables 5 and 6 show the median and best binding energies of distribution results, for the set (test) of 75 docking instances and for the six compared algorithms (gGA, ssGA, DE, and PSO from jMetal; GA and LGA from AutoDock), respectively. These binding scoring values are the energies in kcal/mol associated with the receptor–ligand complex, as previously specified (Section 6.1). Therefore, the lower the binding energy the better the result. Table 8 (in Section 7.2) contains the median RMSD (Root Mean Square Deviation) values associated with each energy binding value obtained. The RMSD is a measure of similarity between the experimental ligand position in the macromolecule and the computed position of the docking ligand. RMSD values are also shown in this study to allow further biological comparisons for the experts in the domain of molecular design. In Table 5, the first observation that can be made is that DE (jMetal) shows the best median energy (marked with gray background) for the highest number of molecular instances. Specifically, DE obtained the best median conformation on 67 out of 75 receptor-ligand complexes, followed by the LGA (AutoDock) with 8 best medians (for molecules 1B6L, 1BDL, 1HEF, 1HIV, 1HPO, 1K6C, 1Z1H, and 1Z1R). In terms of the best binding energy, the DE algorithm also shows the highest number of best solutions (Table 6), achieving the lower binding energy in 37 out of 75 molecules, followed by the LGA with 35. In this regard, the remaining algorithms of jMetal showed the best binding energy in one molecule, obtaining the best solutions: gGA for 1BWB, ssGA for 1Z1H, and PSO for 7HVP. In order to provide these results with statistical confidence (in this study ˛= 0.05), a series of non-parametric statistical tests have been applied, as in several cases the distributions of results did not follow the conditions of normality and homoskedasticity [57]. Therefore, analyses and comparisons focus on the entire distribution of binding energies, although they pay particular attention to the Median values (Table 5), for the 75 tackled molecular instances. In particular, Friedman’s ranking and Holm’s post-hoc multicompare tests [57] have been applied in order to know which algorithms are statistically worse than the control one (the algorithm with the best ranking). In this regard, as shown in Table 7, DE (jMetal) reaches the best ranking value (Friedman) with 1.10, followed by LGA, gGA, ssGA, GA, and PSO. Therefore, DE is Table 3 gGA, ssGA, DE and PSO parameter tuning: Friedman test of all parameter combinations, for each algorithm and for each selected molecular instance (confidence level ˛= 0.05). Each test value in the tables has been calculated from the median distribution out of 30 independent runs. gGA m 5 20 60 100 160 c 510.95 11.27 8.54 6.36 7.22 20 12.54 12.56 14.45 13.09 10.72 60 15.22 13.36 16.09 13.22 15.09 100 14.77 15.04 15.09 16.36 15.36 160 14.04 12.86 12.54 16.13 14.95 According to 2with 24 degrees of freedom: 37.73 With Friedman’s p-value = 3.68E−2 (Rejects H0) DE  0.1 0.3 0.5 0.7 0.9 Cr 0.1 5.90 7.99 9.36 12.63 13.01 0.3 4.54 8.27 13.27 15.90 19.31 0.5 12.09 5.54 12.72 19.00 22.99 0.7 14.72 6.04 7.04 20.90 24.27 0.9 19.09 7.27 2.59 16.63 23.81 According to 2with 24 degrees of freedom: 204.76 With Friedman’s p-value = 1.31E−10 (Rejects H0) ssGA m 520 60 100 160 c 57.99 8.36 12.86 13.68 13.50 20 13.63 17.77 14.27 13.09 15.27 60 13.63 11.95 10.54 12.72 11.95 100 10.95 10.86 14.63 16.13 14.36 160 8.63 16.27 14.68 12.59 14.59 According to 2with 24 degrees of freedom: 29.30 With Friedman’s p-value = 2.01E−1 (Rejects H0) PSO ϕ2 2.05 2.10 2.15 2.20 2.25 ϕ1 2.05 2.40 7.63 5.77 13.63 15.54 2.10 7.09 7.63 12.45 14.81 14.27 2.15 8.13 12.81 12.18 14.18 16.18 2.20 17.36 12.72 16.72 18.72 16.36 2.25 10.04 13.72 18.36 16.27 19.40 According to 2with 24 degrees of freedom: 96.09 With Friedman’s p-value = 2.23E−10 (Rejects H0) Table 4 Parameter settings of compared algorithms. Framework Algorithm Parameter Value jMetal gGA Population size 150 Selection Binary tournament Crossover SBX, pc= 0.9, c= 5 Mutation Polynomial, pm= 1/L,m= 100 ssGA Population size 150 Selection Binary tournament Crossover SBX, pc= 0.9, c= 5 Mutation Polynomial, pm= 1/L,m= 5 DE Population size 150 Differential crossover CR = 0.9 Scaling factor F= 0.5 PSO Swarm size 150 Acceleration coefficients C1= 2.05, C2= 2.05 AutoDock 4.2 GA Population size 150 Crossover 2-point, pc= 0.8 Mutation Cauchy, pm= 0.02, ˛= 0.0, ˇ= 1.0 LGA Population size 150 Crossover 2-point, pc= 0.8 Mutation Cauchy, pm= 0.02, ˛= 0.0, ˇ= 1.0 LS scheme Solis–Wets (300 iterations)