A New Multi-Objective Approach for Molecular Docking Based on RMSD and Binding Energy

A New Multi-Objective Approach for Molecular Docking Based on RMSD and Binding Energy Esteban L´opez-Camacho, Mar´ıa Jes´us Garc´ıa Godoy, Jos´e Garc´ıa-Nieto, Antonio J. Nebro, and Jos´e F. Aldana-Montes Khaos Research Group Departament of Computer Sciences, University of M´alaga, ETSI Inform´atica, Campus de Teatinos, M´alaga - 29071, Spain {esteban,mjgarciag,jnieto,antonio,jfam}@lcc.uma.es Abstract. Ligand-protein docking is an optimization problem based on predicting the position of a ligand with the lowest binding energy in the active site of the receptor. Molecular docking problems are traditionally tackled with single-objective, as well as with multi-objective approaches, to minimize the binding energy. In this paper, we propose a novel multi-objective formulation that considers: the Root Mean Square Deviation (RMSD) difference in the coordinates of ligands and the binding (intermolecular) energy, as two objectives to evaluate the quality of the ligand-protein interactions. To determine the kind of Pareto front approximations that can be obtained, we have selected a set of representative multi-objective algorithms such as NSGA-II, SMPSO, GDE3, and MOEA/D. Their performances have been assessed by applying two main quality indicators intended to measure convergence and diversity of the fronts. In addition, a comparison with LGA, a reference single-objective evolutionary algorithm for molecular docking (AutoDock) is carried out. In general, SMPSO shows the best overall results in terms of energy and and RMSD (value lower than 2 ˚ A for successful docking results). This new multi-objective approach shows an improvement over the ligand-protein docking predictions that could be promising in in silico docking studies to select new anticancer compounds for therapeutic targets that are multidrug resistant. Keywords: Molecular Docking, Multi-Objective Optimization, Nature Inspired Metaheuristics, Algorithm Comparison 1 Introduction Ligand-protein docking is an optimization problem which aims at predicting the position of a small molecule (ligand) to a receptor (macromolecule) with the goal of finding the ligand position to the receptor with a minimum binding energy. Molecular docking problem has been tackled with single-objective algorithms, to minimize the binding energy [11], as well as with multi-objective approaches, to minimize the intermolecular energy Einter (energy interaction between ligand and the target) and the intramolecular energy Eintra (the internal energy compound) [4]. 2 E. L´opez-Camacho et al. In this regard, a number of studies based on the application of multi-objective algorithms to the ligand-protein docking have been proposed. A first attempt was carried out in 2006 by Oduguwa et al. [15], in which three evolutionary multiobjective algorithms (NSGA-II, PAES, and SPEA) were applied to evaluate three objectives such as the Einter , Eintra and shape complementarities on three molecular complexes. Grosdidier et al. [5] proposed a new hybrid evolutionary algorithm called EADock that optimizes two different energy score functions that evaluate the Einter , Eintra and the solvation free energy. In 2008, Janson et al. [7] designed a parallel multi-objective algorithm using AutoDock 3.05 energy function, called ClustMPSO, minimizing as objectives the Einter and Eintra when dealing with six molecular complexes. In the same year, Boisson et al. [1] implemented a parallel evolutionary bi-objective model based on optimizing two objectives: the sum of Einter and Eintra and a surface term for the docking of six instances. Sandoval-Perez et al. [16] used the implementation of NSGA-II provided by the jMetal framework to optimize bound and non-bound energy terms of ligand/receptor as objectives applied to four docking instances. Gu et al. [6] developed a new multi-objective approach based on optimizing the solutions generated by an aggregated scoring function that includes terms from force-field, empirical and knowledge-based scoring functions. In all these previous publications, a series of different multi-objectives formulations have been proposed that focus on energy scoring function. However, they do not consider guiding the search with a new objective when the co-crystallized ligand is known, which could complement the traditional energy function. With this motivation, we propose in this work a novel multi-objective approach consisting minimizing: (1) the binding energy (the unbound and bound energy terms of the ligand/receptor complex), and (2) the Root-Mean-Square-Deviation (RMSD) score, when the co-crystallized ligand pose is known. These two main objectives have been used to evaluate the quality of the ligand-protein interactions. With this aim, we compare and analyze the performance of four multi-objective metaheuristics when solving 11 flexible ligand-receptor docking complexes taken from the AutoDock 4.2 benchmark [12]. This dataset includes flexible ligands with different sizes and flexible side-chains of HIV-protease receptors for more realistic results. The algorithms used in this study are: Nondominated Sorting Genetic Algorithm II (NSGA-II) [2], Speed Modulation Multi-Objective Particle Swarm Optimization (SMPSO) [13], Third Evolution Step of Generalized Differential Evolution (GDE3) [8], and Multi-Objective Evolutionary Algorithm Based on Decomposition (MOEA/D) [18]. These algorithms constitute a varied set of evolutionary and difference-vector multi-objective techniques representative of the state of the art, performing different learning procedures and inducing different behaviors in terms of convergence and diversity. This paper is organized as follows: Section 2 describes the molecular docking problem from a multi-objective formulation. Studied algorithms are briefly described in Section 3. Section 4 reports the experimentation methodology and Section 5 analyzes the results obtained. Finally, Section 6 contains concluding remarks and future lines of research. A new MO Approach for Molecular Docking 3 2 The problem: Multi-Objective Docking A multi-objective optimization problem is characterized by two spaces: the decision and the objective spaces. The former involves all the possible feasible solutions, and the latter includes their corresponding objective values. Decision space. The main objective in the molecular docking problem is to find an optimized conformation between the ligand ( L ) and the receptor ( R ) that results in the lowest binding energy. The ligand-receptor interaction is evaluated by an energy function calculated through three components representing degrees of freedom: (1) the translation of the ligand molecule, involving the three axis values ( x, y, z ) in cartesian coordinate space; (2) the ligand orientation, modeled as a four variables quaternion including the angle slope ( θ ); and (3) the flexibilities, represented by the free rotation of torsion (dihedral angles) of the ligand and sidechains of the receptor. Each problem solution for AutoDock and jMetal (the tools we have used) is encoded by a real-value vector of 7 + n 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 n values are the ligand torsion dihedral angles. Furthermore, in order to allow a rapid evaluation of the energy conformations, a grid-based methodology is implemented. The energy interaction is calculated and assigned to each grid point and is evaluated to obtain the energy of a given ligand pose [12]. Objective space. Our bi-objective formulation consists of: the Einter and the RMSD score. The Einter is the energy function as used in Autodock, that is calculated as follows: Einter =QR−L bound +QR−L unbound (1) QR−L bound and QR−L bound are the states of bound and unbound of the ligand-receptor complex, respectively. Q=Wvdw X i,j (Aij r12 ij − Bij r6 ij ) + Whbond X i,j E(t) Cij r12 ij − Dij r10 ij !+ +Welec X i,j qiqj ε(rij)rij +Wsol X i,j (SiVj+SjVi)e(−r2 ij /2σ2) (2) Each pair of energetic evaluation terms includes evaluations ( Q ) of dispersion/repulsion ( vdw ), hydrogen bonds ( hbond ), electrostatics ( elec ) and desolvation ( sol ). Weights Wvdw , Whbond , Wconf , Welec , and Wsol of Equation 2 are constants for Van der Waals, hydrogen bonds, torsional forces, electrostatic interactions and desolvation, respectively. rij represents the interatomic distance, Aij and Bij in the first term are Lennard-Jones parameters taken from the Amber force field. 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 uses a 4 E. L´opez-Camacho et al. 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. An extended explanation of all these variables can be found in [12]. The RMSD is a measure of similarity between the real ligand position in the receptor and the computed position of the docking ligand, that takes into account symmetry, partial symmetry (e.g. symmetry within a rotatable branch) and near-symmetry in a simple heuristic way. Ideally, the lower RMSD score the better solution is. A ligand-receptor docking solution with a RMSD score below 2 ˚ A is considered as a solution with high docking accuracy. It is worth noting that other docking solutions can be returned with higher RMSD scores and low values of Einter , indicating that other possible interaction ligand sites should be considered. The RMSD score for two identical structures a and b is defined as follows: RMSDab =max(RMSD0 ab, RMSD0 ba), with RMSD0 ab =s1 NX i min jrij 2 (3) The sum is over all N heavy atoms in structure a , the minimum is over all atoms in structure awith the same element type as atom iin structure b. 3 Algorithms We have included in our study four algorithms which are representative of the state-of-the-art in the multi-objective optimization field. A brief description of each one of them is given next: NSGA-II: NSGA-II [2] is a generational genetic algorithm, which uses the typical genetic operators (selection, crossover and mutation) to obtain new individuals from the original population. To promote convergence, a non-dominated sorting procedure based on Pareto ranking is used, while the crowding distance density estimator is applied to foster the diversity of the set of found solutions. GDE3: The Generalized Differential Evolution (GDE) algorithm [8] is based on NSGA-II, but the genetic mutation and selection operators are replaced by their differential evolution counterparts. Furthermore, GDE3 modifies the crowding distance of NSGA-II as well to generate a better distributed set of solutions. SMPSO: SMPSO [13] is a multi-objective particle swarm optimization algorithm. Its main feature is the limitation of the particle speed to allow new effective particle positions to be produced when the speed becomes too high. SMPSO uses the polynomial mutation as the turbulence factor and an external archive that stores the non-dominated solutions found during the search. A new MO Approach for Molecular Docking 5 MOEA/D: MOEA/D [18] has become the typical representative decompositionbased multi-objective algorithm, where a multi-objective problem is decomposed into a set of single-objective subproblems that then optimized simultaneously. In this study we have used the variant MOEA/D-DE [9], which applies differential evolution as variation operators. This algorithm also applies a polynomial mutation operator to improve its search capability. In short, we have selected the most widely used algorithm in the field (NSGAII), a solver based in differential evolution (GDE3), a PSO (SMPSO) and an algorithm based on decomposition (MOEA/D). 4 Experimentation In this section, we include the selected benchmark problems, the experimentation methodology we have followed, and the parameter settings of the algorithms. 4.1 Benchmark Problems In this study, we have selected a benchmark composed of 11 complexes having receptor and ligand flexibility. The selection of these complexes has been motivated as they are actually difficult docking problems containing a wide range of ligand sizes (from small to large inhibitors). The receptors of these complexes have a tunnel-shaped active site that wraps around a peptidomimetic inhibitor [12]. The receptor is a dimer whose subunits are bridged by an arginine-aspartate salt bridge at the end of the tunnel. The docking studies performed with these instances in [12] to test the energy function of AutoDock 4.2 demonstrated that the most difficult problems are those which involve smaller ligands. This is due to the flexibility added to the receptor side-chains (ARG-8) that increases the space of ligand interaction. These instances have been taken from the PDB database 1 and they have been properly prepared for the docking simulations. Table 1 summarizes the set of problems selected showing the PDB accession code, the X-ray crystal structures names and the structure resolution ( ˚ A). For all instances, the torsional degrees of freedom (flexibility) for ligands and macromolecules are 10 and 6, respectively, selecting those torsions that allow the fewest number of atoms to move around the ligand core. Therefore, the total number of solution variables ( n ) is 23 (3 for translation, 4 for rotation quaternion, and 16 for torsional degrees). 4.2 Methodology For this work, we have carried out a thorough experimentation consisting in performing 30 independent runs for each combination of algorithm and molecular instance. From these executions, we have calculated the median and interquartile range (IQR) as measures of central tendency and statistical dispersion, 1In URL: http://www.rcsb.org/pdb/home/home.do 6 E. L´opez-Camacho et al. Table 1. The accession codes, the X-ray crystal structure and resolution taken from PDB database are presented. PDB Code Protein-ligand complexes Resolution (˚ A) 1AJV HIV-1 protease/AHA006 2.00 1AJX HIV-1 protease/AHA001 2.00 1BV9 HIV-1 protease/α-D-glucose 2.20 1D4K HIV-1 protease/Macrocyclic peptidomimetic inhibitor 8 1.85 1G2K HIV-1 protease/AHA047 1.95 1HIV HIV-1 protease/U75875 2.00 1HPX HIV-1 protease/KNI-272 2.00 1HTF HIV-1 protease/GR126045 2.20 1HTG HIV-1 protease/GR137615 2.00 1HVH HIV-1 protease/Q8261 1.80 2UPJ HIV-1 protease/U100313 3.00 respectively. We have considered two quality indicators to assess the algorithm performance: Hypervolume ( IHV ) and Unary Additive Epsilon Indicator ( I+ ) [3]. The first indicator takes into account both convergence and diversity, whereas the second one ( I+ ) gives a measure of the convergence degree of the obtained Pareto front approximations. In this sense, it is worth noting that we are dealing with a real-world optimization problem, and therefore the true Pareto fronts to calculate these two metrics are not known. To cope with this issue, we have generated a reference Pareto front for each instance by combining all the non-dominated solutions computed in all the executions of all the algorithms. As mentioned, we have used the implementation of the four algorithms studied provided in the jMetalCpp framework [10] in combination with AutoDock 4.2 to evaluate the new generated solutions. To cope with the high computational requirements needed by carry out our experiments, we have used the Condor 2 system, a middleware platform managing close to 400 cores that acts as a distributed task scheduler (each task dealing with one independent run). 4.3 Parameter Setup The selected algorithms have been configured with a population size of 150 individuals (particles in the case of SMPSO). The stopping condition has been set to compute a number of 1,500,000 function evaluations. These values were chosen because they are the default settings used by AutoDock and they have been used in other studies [14]. Each algorithm has been configured using the parameter setup recommended in the research study where it was proposed, and these parameters are used as default in the jMetal framework. In particular, SBX crossover and polynomial mutation are the variation operators used in NSGA-II. The distribution indexes for both operators are ηc = 20 for crossover, and ηm = 20 for mutation. The 2In URL: http://research.cs.wisc.edu/htcondor/ A new MO Approach for Molecular Docking 7 Table 2. Median and interquartile range of IHV for each algorithm and instance. Best and second best median results have dark and light gray backgrounds, respectively. NSGAII SMPSO GDE3 MOEAD 1AJV 0.00e+ 000.0e+00 3.51e−014.0e−02 0.00e+ 000.0e+00 0.00e+ 002.9e−01 1AJX 0.00e+ 000.0e+00 5.52e−012.0e−02 0.00e+ 000.0e+00 7.47e−036.8e−01 1D4K 0.00e+ 000.0e+00 4.93e−011.3e−01 0.00e+ 000.0e+00 0.00e+ 000.0e+00 1G2K 0.00e+ 000.0e+00 3.32e−013.3e−02 0.00e+ 000.0e+00 0.00e+ 004.1e−01 1HIV 0.00e+ 000.0e+00 5.96e−011.3e−01 0.00e+ 000.0e+00 0.00e+ 000.0e+00 1HPX 0.00e+ 000.0e+00 2.04e−011.8e−01 1.27e−016.5e−01 0.00e+ 001.1e−01 1HTF 0.00e+ 000.0e+00 5.26e−021.3e−01 0.00e+ 004.6e−03 2.78e−023.3e−01 1HTG 0.00e+ 000.0e+00 3.51e−025.6e−02 0.00e+ 000.0e+00 0.00e+ 001.9e−01 1HVH 0.00e+ 000.0e+00 7.67e−013.7e−02 0.00e+ 000.0e+00 5.31e−017.7e−01 1VB9 0.00e+ 000.0e+00 7.34e−016.5e−02 0.00e+ 000.0e+00 0.00e+ 001.4e−01 2UPJ 0.00e+ 000.0e+00 5.86e−019.8e−02 0.00e+ 000.0e+00 1.90e−015.8e−01 crossover probability is pc = 0 . 9 and the mutation probability is pm = 1 /n , being n the number of decision variables of the tackled problem. NSGA-II applies binary tournament selection. In the case of GDE3 (variant rand/1/bin), the two DE control parameters µ and Cr take a value of 0.5, whereas in MOEA/D µ is set to 0.5 and Cr is set to 1.0. Both MOEA/D and SMPSO use the polynomial mutation with the same settings applied in NSGA-II. In SMPSO, the acceleration coefficients ϕ1 and ϕ2 are set to 1.5, the inertia weight is W = 0 . 9, and the polynomial mutation is applied to one sixth of the particles in the swarm. 5 Results This section is devoted to presenting and analyzing the results obtained in our study. We start by assessing the performance of the algorithms and then they are compared with the values of a single-objective approach. 5.1 Performance Comparisons We start our analysis by discussing the results yielded by applying the IHV indicator. Let us remind that this indicator is the sum of the contributed volume of each point of a front in respect to a reference point, and the higher the convergence and diversity degrees of a front, the higher its IHV value. Table 2 shows the median and interquartile range of the computed solutions for IHV quality indicators for the set of 11 docking instances and the four algorithms being compared. According to these results, SMPSO achieves the best IHV values in all the eleven considered problems and MOEA/D is the second best performing technique. We have to note that many cells have a IHV value equal to zero; this happens when all the points of the produced fronts are beyond the limits of the reference point. This happens in most of the problems in all the algorithms excepting SMPSO, which indicates we are facing a very hard optimization problem. A similar behavior can be observed in Table 3 with regards to I+ , which is a convergence measure. According to these results, SMPSO obtains the best I+ 8 E. L´opez-Camacho et al. Table 3. Median and interquartile range of I+ for each algorithm and instance. Best and second best median results have dark and light gray backgrounds, respectively. NSGAII SMPSO GDE3 MOEAD 1AJV 5.23e+ 001.2e+00 5.60e−019.8e−02 5.00e+ 001.0e+00 3.87e+ 004.4e+00 1AJX 3.43e+ 002.4e+00 2.61e−017.2e−02 1.49e+ 003.3e−01 1.01e+ 002.0e+00 1D4K 8.06e+ 002.7e+00 4.56e−011.4e−01 8.56e+ 005.7e−01 4.65e+ 002.8e+00 1G2K 4.28e+ 001.4e+00 5.71e−011.2e−01 3.93e+ 001.3e+00 2.69e+ 003.6e+00 1HIV 5.12e+ 001.2e+00 2.63e−012.1e−01 4.69e+ 001.4e+00 4.07e+ 001.6e+00 1HPX 1.42e+ 013.6e+00 6.32e−012.8e−01 6.71e−011.1e+01 1.03e+ 011.3e+01 1HTF 1.76e+ 005.5e−01 9.30e−013.0e−01 1.13e+ 008.0e−01 7.94e−019.2e−01 1HTG 7.48e+ 007.1e−01 9.63e−016.6e−02 6.82e+ 008.7e−01 5.03e+ 006.4e+00 1HVH 5.94e+ 001.5e+00 1.34e−012.7e−02 4.93e+ 001.7e+00 4.16e−012.1e+00 1VB9 8.59e+ 002.4e+00 1.33e−015.6e−02 7.85e+ 001.3e+00 7.04e+ 004.9e+00 2UPJ 3.42e+ 002.4e+00 3.03e−016.6e−02 3.56e+ 001.1e+00 7.64e−012.7e+00 Table 4. Average Friedman’s rankings with Holm’s Adjusted p -values ( α = 0 . 05) of compared algorithms (SMPSO, GDE3, MOEA/D, and NSGA-II) for the test set of 11 docking instances. Symbol * indicates the control algorithm and column at right contains the overall ranking of positions with regards to IHV and I+. Hypervolume (HV) Epsilon (I+) Algorithm F riRank HolmAp Algorithm F riRank HolmAp *SMPSO 1.02 -*SMPSO 1.09 - MOEA/D 2.68 2.24e-03 MOEA/D 2.00 9.87e-02 GDE3 3.09 1.45e-04 GDE3 3.09 2.79e-04 NSGA-II 3.22 5.21e-05 NSGA-II 3.81 7.25e-07 values in ten out of the eleven problems, while MOEA/D gets the best value in one problem (1HTF) and the second best in all of them but 1HPX. In order to provide these results with statistical meaning (in this study α = 0 . 05), non-parametric statistical tests have been applied because in several cases the distributions of results did not follow the conditions of normality and homoscedasticity [17]. Therefore, the analyses and comparisons focus on the entire distribution of each of the two metrics studied. Specifically, we have applied Friedman’s ranking and Holm’s post-hoc multicompare tests [17] to know which algorithms are statistically worse than the control one (with the best ranking). In this regard, as shown in Table 4, SMPSO is the best ranked technique according to IHV (with a value of 1.02), followed by MOEA/D, GDE3, and NSGA-II. Therefore, SMPSO is established as the control algorithm for IHV in the post-hoc Holm test, which is compared with the remaining algorithms. The adjusted p-values ( HolmAp in Table 4) resulting from these comparisons are, for the last three algorithms (MOEA/D, GDE3, and NSGA-II), lower than the confidence level, meaning that SMPSO is statistically better than these algorithms. In the case of I+ , SMPSO is better ranked than the remaining compared algorithms, although without statistical differences in the case of MOEA/D. SMPSO is statistically better than GDE3 and NSGA-II. In summary, SMPSO shows the overall best balance for the two quality indicators, followed by MOEA/D. These results are graphically supported by two examples included in Fig. 1, where the reference fronts obtained for two A new MO Approach for Molecular Docking 9 −30 −29 −28 −27 0.0 0.2 0.4 0.6 1D4K Energy (kcal/mol) RMSD (Å) MOEA/D RF −30 −29 −28 −27 0.0 0.2 0.4 0.6 1D4K Energy (kcal/mol) RMSD (Å) SMPSO RF −27.0 −26.5 −26.0 −25.5 −25.0 −24.5 0.0 0.2 0.4 0.6 0.8 1.0 1HIV Energy (kcal/mol) RMSD (Å) MOEA/D RF −27.0 −26.5 −26.0 −25.5 −25.0 −24.5 0.0 0.2 0.4 0.6 0.8 1.0 1HIV Energy (kcal/mol) RMSD (Å) SMPSO RF Fig. 1. Reference front contributions of docking instances 1D4K and 1HIV. SMPSO and MOEA/D contribute with practically all the solutions of the reference fronts. representative instances 1D4K and 1HIV are plotted. In these graphs, the contributions, in terms of solutions, of each algorithm to the global reference front are plotted with different points and colors. As it is easily observable, SMPSO and MOEA/D contribute with almost all solutions taking part to the reference front. Interestingly, SMPSO converges to the region biased towards the RMSD objective, whereas MOEA/D generate non-dominated solutions in a different region to the ones of SMPSO, thereby giving cue to the energy optimization. We can state that the specific learning procedures induced by SMPSO and MOEA/D lead these algorithms to search in different regions of the problem landscape, hence generating solutions in complementary parts of the reference front. 5.2 Comparison Single Versus Multi-objective After the performance comparison of multi-objective algorithms, we are now interested in knowing how competitive their solutions are against those yielded