A Study of Archiving Strategies in Multi-Objective PSO for Molecular Docking

A Study of Archiving Strategies in Multi-Objective PSO for Molecular Docking Jos´e Garc´ıa-Nieto1, Esteban L´opez-Camacho1, Mar´ıa Jes´us Garc´ıa Godoy1, Antonio J. Nebro1, Juan J. Durillo2, and Jos´e F. Aldana-Montes1 1Khaos Research Group Departament of Computer Sciences, University of M´alaga, ETSI Inform´atica, Campus de Teatinos, M´alaga, Spain jnieto,esteban,mjgarciag,antonio,[email protected] 2Distributed and Parallel Systems Group, University of Innsbruck, Austria [email protected] Abstract. Molecular docking is a complex optimization problem aimed at predicting the position of a ligand molecule in the active site of a receptor with the lowest binding energy. This problem can be formulated as a bi-objective optimization problem by minimizing the binding energy and the Root Mean Square Deviation (RMSD) difference in the coordinates of ligands. In this context, the SMPSO multi-objective swarm-intelligence algorithm has shown a remarkable performance. SMPSO is characterized by having an external archive used to store the non-dominated solutions and also as the basis of the leader selection strategy. In this paper, we analyze several SMPSO variants based on different archiving strategies in the scope of a benchmark of molecular docking instances. Our study reveals that the SMPSOhv, which uses an hypervolume contribution based archive, shows the overall best performance. Key words: Multi-Objective Optimization, Particle Swarm Optimization, Molecular Docking, Archiving Strategies, Algorithm Comparison 1 Introduction Molecular docking is a complex optimization problem found in biology, which consists in predicting the position of a small molecule (ligand) in the active site of a receptor (macromolecule) that registers the minimum binding energy. Molecular docking is traditionally faced by means of metaheuristics [8] [4] as a continuous optimization problem, since it requires to adjust position variables corresponding to coordinates of translation and torsion movements of molecules. In the last decade, a number of studies have centered on the application of singleand multi-objective metaheuristics [8] [4] [6] [15] [5] to the molecular docking problem, showing successful results for a number of molecular compounds. In these previous works, different objective formulations were proposed that focused on energy scoring functions. Recently, a new multi-objective approach has been proposed [9] in which two different objectives are to be minimized: the 2 Jos´e Garc´ıa-Nieto et al. binding energy (the unbound and bound energy terms of the ligand/receptor complex) and the Root-Mean-Square-Deviation (RMSD) score. The latter objective leads the algorithms to guide the search when the co-crystallized ligand is known, which complements the traditional energy function. Among these optimization techniques, a multi-objective swarm-intelligence approach, namely SMPSO [11], has emerged as one of the most prominent optimizers for molecular docking [4] [9]. This technique performs a limitation mechanism of particle’s velocity to avoid the movement of particles in search regions out of the problem ranges. SMPSO uses an external archive to store non-dominated solutions according to the crowding distance [2]. This archive is also used in the leader selection mechanism. Here, our motivation is to go one step beyond by evaluating, in the scope of a benchmark of molecular instances, new versions of SMPSO using different archiving strategies (hypervolume, cosine distance, and aggregation) and, consequently, different strategies for the selection of the leaders. With this aim, we compare and analyze the proposed versions of SMPSO when solving 11 flexible ligand-receptor docking complexes taken from the AutoDock 4.2 benchmark [10]. This dataset includes flexible ligands with different sizes and flexible side-chains of HIV-protease receptors.The performance of the algorithms has been assessed by applying two main quality indicators intended to measure convergence and diversity of the computed Pareto front approximations. The remainder of this article is organized as follows: Section 2 describes the molecular docking problem from a multi-objective formulation. Studied algorithms are described in Section 3. Section 4 reports the experimentation methodology and Section 5 analyzes the obtained results. Finally, Section 6 reports conclusions and future lines of research. 2 Molecular Docking From a biological point of view, 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 a minimum binding energy. The interaction between Land Rcan be described by an energy function calculated from 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. - Solution Encoding: Each problem solution is then encoded by a realvalue vector of 7 + nvariables (as illustrated in Fig. 1), in which the first three values correspond to the ligand translation, the next four values correspond to the ligand and/or receptor orientation, and the remaining nvalues are the ligand torsion dihedral angles. The range of translation variables (x, y, z) is [0 · · · 120], which has been delimited between the limits of the coordinates of a grid space previously set for On the Evaluation of Archiving Strategies in SMPSO 3 Fig. 1. Solution encoding. The first three values (translation) are the coordinates of the center of rotation of the ligand. The next four values correspond to quaternion and (θ). The rest of the values hold the torsion angles in degrees. each problem. Orientation (quaternion) and torsion variables are measured in radians and encoded in the range of [−π, π]. - Fitness Functions: the bi-objective formulation used here consists of: the Ebinding and the RMSD score. The Ebinding is the energy function as used in Autodock, which is calculated as follows: Ebinding =QR−L bound +QR−L unbound (1) 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) QR−L bound and QR−L bound are the states of bound and unbound of the ligand-receptor complex, respectively. 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. An extended explanation of all these variables can be found in [10]. The RMSD is a measure of similarity between the real ligand position in the receptor and the computed position of the docking ligand. The lower RMSD score the better the solution is. A ligand-receptor docking solution with an RMSD score below 2˚ A is considered as a solution with high docking accuracy. 4 Jos´e Garc´ıa-Nieto et al. The RMSD score for two identical structures aand bis defined as follows: RMSDab =max(RMSD 0 ab, RMSD 0 ba), with RMSD 0 ab =s1 NX i min jrij 2 (3) The sum is over all Nheavy atoms in structure a, the minimum is over all atoms in structure awith the same element type as atom iin structure b. 3 Algorithms In this section, we describe the SMPSO variants we are going to study. We start with the original algorithm and then we give details of the considered variants. SMPSO is a Multi-Objective Particle Swarm Optimization (MOPSO) characterized by two features: a velocity constraint mechanism and an external bounded archive to store the non-dominated solutions found during the search [11]. A perturbation, implemented as a mutation operator, is also incorporated. Its pseudocode is included in Algorithm 1. The archive contains the current Pareto front approximation found by the algorithm, and it applies the crowding distance density estimator [2] to decide which particle to remove when it is full. The archive is also used in the leader strategy selection, consisting on binary tournament based on randomly selecting two solutions from it and taking the one with the highest crowding distance value (i.e., the one located in less crowded region of the front composed by all archived solutions). The local best position of a particle iis obtained by applying a dominance test with the rest of particles in the swarm, in such a way that the current best particle (which initially is particle i) is updated when it is dominated by another one. In [12], a study of different leader selection mechanisms on SMPSO was conducted. In that work, the most salient variant consisted in replacing the crowding distance by the degree of contribution of the solutions in the external archive according to the hypervolume indicator [18]. This way, the leader selection is based on a binary tournament that chooses the particle having the largest hypervolume contribution value. This version was named as SMPSOhv and it is the second selected algorithm to be compared in our study. We introduce in this paper a new variant of SMPSO. The cosine similarity is a measure of similarity between two vectors that measures the cosine of the angle between them. This way, two vectors in the same direction have a cosine similarity value equals to zero, while two perpendicular vectors have a cosine similarity value of 1. As all the solutions in an external archive are non-dominated, we can define a density estimator by fixing a reference point and computing the cosine similarity among the vectors conformed by the archive solutions with regards to that reference point. The studied problem in this paper has two objectives, so we can sort the solutions in the archive by the first objective and compute, for each solution, a density value by summing up the cosine similarity of each point to their previous and next points; extreme points have a similarity distance equals to 0. This way, points having a largest cosine density value are in On the Evaluation of Archiving Strategies in SMPSO 5 Algorithm 1 Pseudocode of SMPSO 1: initializeSwarm() 2: initializeLeadersArchive() 3: generation = 0 4: while generation <maxGenerations do 5: computeSpeed() 6: updatePosition() 7: mutation() // perturbation 8: evaluation() 9: updateLeadersArchive() 10: updateParticlesMemory() 11: generation ++ 12: end while 13: returnLeadersArchive() the most densely populated region. The resulting algorithm is called SMPSOC. An important issue in this technique is to select the proper reference point. Our previous study [9] indicated that the fronts have a convex shape, so we choose an approximation to the nadir point by taking the highest objective values of the solutions in the archive. The fourth SMPSO version in our study, also presented in this paper for the first time, is an archive-less approach and it is called SMPSOD. To leave out the archive, we take the strategy of designing an aggregative version of SMPSO inspired by MOEA/D [17], where a multi-objective problem can be decomposed into a set of single-objective problems that can be optimized at the same time. This way, a set of evenly spread weight vectors λ1, λ2, . . . , λNare defined, being Nthe size of the swarm. Then, each particle ihas associated the vector λiand a neighborhood defined as a set of its several closest weight vectors in λ1, λ2, . . . , λN. The scalarizing strategy follows the Tchebycheff scheme. The strategy for getting the local best of a particle iis the same procedure as used by MOEA/D to update a neighborhood. The leader updating strategy consists in finding the best solution in the neighborhood by considering the scalar values of the particles taking into account their weight vectors. SMPSO was inspired in the OMOPSO algorithm proposed by Reyes and Coello in [14], so we decided to also include it in our comparisons as a reference multi-objective particle swarm optimizer in the state of the art. In summary, in our study we include OMOPSO and four SMPSO variants with different archiving strategies: crowding distance based (original SMPSO), hypervolume contribution based (SMPSOhv), cosine distance based (SMPSOC), and archive-less (SMPSOD), being the last two ones proposed in this paper. 4 Experimentation For the experiments, we have considered a benchmark of 11 molecular instances with receptor and ligand flexibility. These complexes are actually difficult dock- 6 Jos´e Garc´ıa-Nieto 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 ing problems containing a wide range of ligand sizes (from small to large inhibitors). The docking studies performed with these instances in [10] to test the energy function of AutoDock 4.2 demonstrated that the most difficult problems are those involving smaller ligands. This is due to the flexibility added to the receptor side-chains (ARG-8) that increases the space of ligand interactions. These instances have been taken from the PDB database3. 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 for ligands and receptors are 10 and 6, respectively, selecting those torsions that allow the fewest number of atoms to move around the ligand core. Therefore, the solution vector contains: 3 variables for translation, 4 variables for rotation quaternion, and 16 variables for torsional degrees, summing up a total number (n) of 23 variables. 4.1 Methodology The followed methodology consists in running each combination of algorithm and molecular instance 30 independent times. From these executions, we have calculated the median and interquartile range (IQR) as measures of central tendency and statistical dispersion, respectively. We have considered two quality indicators to assess the algorithm performance: Hypervolume (IHV ) [18] and Unary Additive Epsilon Indicator (I+) [19]. The former takes into account both convergence and diversity, whereas the later gives a measure of the convergence degree of the obtained Pareto front approximations. It is worth noting that we are dealing with a real-world optimization problem, and therefore the 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. We have used the implementation of the five studied algorithms provided in the jMetalCpp framework [7], in combination with AutoDock 4.2 to evaluate the new generated solutions. To cope with the high computational requirements 3In URL: http://www.rcsb.org/pdb/home/home.do On the Evaluation of Archiving Strategies in SMPSO 7 Table 2. Parameter settings. Common parameters Swarm size 150 Particles Iterations 10,000 SMPSO [3] & SMPSOhv & SMPSOD & SMPSOC Archive Size 100 C1,C21.5 w0.9 Mutation polynomial mutation Mutation probability 1.66 Mutation distribution index ηm20 Selection method Rounds OMOPSO [1] Archive size 100 C1,C2rand(1.5,2.0) w rand(0.1,0.5) Mutation uniform + non-uniform + no mutation Mutation probability Each mutation is applied to 1/3 of the swarm needed to carry out all the experiments, we have used the Condor4system, a middleware platform acting a distributed task scheduler of up to 400 cores. The parameter settings are summarized in Table 2. We set a common subset of parameters which are the same for all the evaluated algorithms. The size of the swarm is 150 and the stopping condition is reached when 1,500,000 function evaluations are performed. These values were chosen as they are the default settings in AutoDock and they have been used in previous studies [13]. The archive size, when applicable, is set to 100. All SMPSO versions use the polynomial mutation with distribution index ηm= 20, which is applied to one sixth of the particles in the swarm. The acceleration coefficients C1and C2are set to 1.5 and the inertia weight is w= 0.9. With these parameters setting, our approach has been to use common settings in order to make a fair comparison, keeping the rest of the parameters of SMPSO and OMOPSO according to the papers where they were originally described. 5 Results and Analysis A first analysis in our experimentation corresponds to the of results in terms of the hypervolume indicator IHV . This indicator computes the sum of the contributed volume of each point in the Pareto front (non-dominated solutions) with regards to a reference point. Therefore, the higher the convergence and diversity degree of a front, the higher (better) the resulting IHV value is. Table 3 shows the median and interquartile range of the computed distributions (out of 30 independent runs) of IHV , for the set of 11 docking instances and for the five compared algorithms. As we can observe, SMPSOhv obtains the best median values of IHV for all the molecular instances and SMPSOD is the second best performing technique. We have to mention that some results of OMOPSO have a IHV equal to zero. This happens when all the points of the produced fronts are dominated by the reference point. In contrast, all the 4In URL: http://research.cs.wisc.edu/htcondor/ 8 Jos´e Garc´ıa-Nieto et al. Table 3. Median and interquartile range of IHV for each algorithm and instance. Best and second best median results have dark and light gray backgrounds, respectively. SMPSO SMPSOhv SMPSOD SMPSOC OMOPSO 1AJV 3.65e−015.1e−02 4.33e−014.0e−02 3.63e−014.6e−02 3.55e−014.8e−02 0.00e+ 000.0e+00 1AJX 4.31e−012.5e−02 5.06e−012.7e−02 4.74e−013.7e−02 4.43e−013.6e−02 0.00e+ 000.0e+00 1D4K 6.67e−018.1e−02 8.48e−011.1e−01 7.11e−019.4e−02 7.35e−019.2e−02 0.00e+ 000.0e+00 1G2K 3.84e−015.3e−02 4.58e−015.9e−02 3.82e−014.1e−02 3.52e−015.2e−02 0.00e+ 000.0e+00 1HIV 4.86e−012.0e−01 6.74e−012.9e−02 5.87e−017.1e−02 4.66e−012.4e−01 0.00e+ 000.0e+00 1HPX 3.60e−011.8e−01 6.30e−019.7e−02 4.77e−011.0e−01 4.63e−011.4e−01 0.00e+ 000.0e+00 1HTF 2.61e−013.3e−01 4.17e−012.4e−01 3.96e−017.9e−02 2.77e−013.1e−01 0.00e+ 000.0e+00 1HTG 8.33e−021.3e−01 1.46e−019.6e−02 1.03e−018.2e−02 7.13e−021.3e−01 0.00e+ 000.0e+00 1HVH 7.78e−014.7e−02 8.69e−019.3e−03 7.70e−012.4e−02 7.85e−012.9e−02 0.00e+ 000.0e+00 1VB9 4.10e−011.2e−01 5.09e−015.6e−02 4.12e−011.1e−01 4.38e−019.1e−02 0.00e+ 000.0e+00 2UPJ 5.82e−019.6e−02 6.96e−015.1e−02 6.27e−017.4e−02 6.20e−016.8e−02 1.99e−016.4e−01 Table 4. Median and interquartile range of I+for each algorithm and instance. Best and second best median results have dark and light gray backgrounds, respectively. SMPSO SMPSOhv SMPSOD SMPSOC OMOPSO 1AJV 5.12e−011.0e−01 3.94e−016.7e−02 5.35e−011.0e−01 5.46e−011.0e−01 5.31e+ 002.0e+00 1AJX 2.31e−011.1e−01 1.32e−014.3e−02 1.94e−016.1e−02 2.57e−019.4e−02 2.54e+ 003.2e+00 1D4K 2.06e−018.6e−02 4.41e−021.2e−01 1.54e−017.3e−02 1.57e−018.1e−02 8.81e+ 004.1e+00 1G2K 4.29e−011.7e−01 2.81e−012.0e−01 4.75e−019.7e−02 5.15e−011.1e−01 6.01e+ 002.3e+00 1HIV 3.95e−013.6e−01 9.03e−026.4e−02 2.66e−011.2e−01 4.36e−013.2e−01 4.91e+ 001.1e+00 1HPX 4.25e−012.8e−01 1.30e−019.2e−02 2.95e−011.2e−01 3.17e−011.7e−01 1.13e+ 015.7e+00 1HTF 6.60e−011.5e+00 5.46e−013.7e−01 5.64e−011.1e−01 6.85e−014.3e−01 1.49e+ 006.2e−01 1HTG 9.07e−011.4e−01 8.35e−019.3e−02 8.84e−018.8e−02 9.23e−012.1e−01 1.21e+ 017.7e+00 1HVH 1.46e−014.4e−02 6.12e−024.8e−03 1.47e−013.9e−02 1.52e−013.1e−02 5.11e+ 002.4e+00 1VB9 3.34e−012.2e−01 1.96e−017.7e−02 3.44e−011.8e−01 2.97e−011.3e−01 9.31e+ 001.6e+00 2UPJ 2.86e−017.9e−02 1.76e−019.2e−02 2.25e−011.4e−01 2.70e−015.1e−02 7.74e−014.0e+00 SMPSO versions obtained IHV values higher than zero, which indicates that they are all able to produce solutions within the limits of the reference point. In the case of I+, a similar observation can be extracted from Table 4. That is, SMPSOhv shows the best results for all instances, followed by SMPSOD and SMPSO (the lower I+value, the better the result is). For this indicator, SMPSOC obtains a second best median value only for instance 1VB9. These results are assessed with statistical confidence (in this study p-value = 0.05) by focusing on the entire distribution of each of the two studied metrics. In concrete, we have applied Friedman’s ranking and Holm’s post-hoc multicompare tests [16] to know which algorithms are statistically worse than the control one (i.e, the one ranking the best). This way, as shown in Table 5, SMPSOhv is the best ranked variant according to Friedman test for the two indicators (IHV and I+), and it is followed by SMPSOD. Therefore, SMPSOhv is established as the control algorithm in the post-hoc Holm tests, which is compared with the rest of algorithms. The adjusted p-values (HolmAp in Table 5) resulting from these comparisons are, for the remaining variants (SMPSO, SMPSOD, SMPSO, and OMOPSO), lower than the confidence level (0.05), meaning that SMPSOhv is statistically better than these algorithms. SMPSO and SMPSOC obtained similar overall performances, although showing SMPSOC better ranking than SMPSO in terms of IHV . On the Evaluation of Archiving Strategies in SMPSO 9 Table 5. Average Friedman’s rankings with Holm’s Adjusted p-values (0.05) of compared algorithms 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 (IHV ) Epsilon (I+) Overall Algorithm F riRank HolmAp Algorithm F riRank HolmAp Algorithm Rank *SMPSOhv 1.01 -*SMPSOhv 1.00 -SMPSOhv 2 SMPSOD 2.54 2.18e-02 SMPSOD 2.45 3.09e-02 SMPSOD 4 SMPSOC 3.09 3.85e-03 SMPSO 2.99 6.02e-03 SMPSOC 5 SMPSO 3.36 1.36e-03 SMPSOC 3.54 4.79e-04 SMPSO 5 OMOPSO 4.99 1.19e-08 OMOPSO 4.98 1.19e-08 OMOPSO 10 Fig 2 shows the boxplots of the distributions of results concerning the IHV values, for each compared algorithm and molecular instance. In this figure, we can check that SMPSOhv variant obtains the best distributions for all the instances. An interesting observation can be made regarding OMOPSO, whose distributions denote poor results, although it produces outlier solutions with the best indicator values for some instances: 1AJV, 1AJX, 1HPX, and 1HTF. These outliers lead OMOPSO to contribute with many solutions to the reference Pareto fronts. An example of this can be observed in Fig. 3, where the fronts with best IHV values of all compared algorithms are plotted for instance 1AJX. However, the overall results (in boxplots) of OMOPSO indicate that it behaves irregular (non-robust) for all the molecular instances. Following with Fig. 3, another interesting observation lies in the ability of SMPSOD to obtain non-dominated solutions in the region of the reference Pareto front with low energy and high RMSD values (top-left in plots of Fig. 3). In contrast with the other compared algorithms, SMPSOD is able to properly cover this area, as well as other areas with low RMSD. Therefore, as suggested in our previous study [9], a hybrid implementation of SMPSO using an aggregative (archive-less) strategy as done in MOEA/D, would cover the reference front with non-dominated solutions in the two objective ends. This assumption is now tested with SMPSOD in this study. In summary, SMPSOhv shows the overall best behaviour followed by SMPSOD. Intuitively, the former obtains the best IHV as it performs a leader selection method of non-dominated solutions (from the external archive) with largest hypervolume contributions. That is, the particles in the swarm are guided by leaders with large hypervolume contributions, which would enable SMPSOhv to obtain, not only high values of IHV , but also accurate results in terms of I+. 6 Conclusions In this paper, we analyze new variants of SMPSO, a multi-objective swarm optimization technique, based on different archiving strategies in the scope of a benchmarking set of molecular docking instances. The problem is formulated as a bi-objective optimization problem, by minimizing the binding energy and the Root Mean Square Deviation (RMSD) difference in the coordinates of ligands.