Bio-inspired optimization for the molecular docking problem: State of the art, recent results and perspectives

Bio-inspired optimization for the molecular docking problem: State of the art, recent results and perspectives María Jesús García-Godoy a, Esteban López-Camacho a, José García-Nieto a, Javier Del Ser b,c,d,∗, Antonio J. Nebro a, José F. Aldana-Montes a a University of Málaga, 29071 Málaga, Spain b TECNALIA, 48160 Derio, Spain c University of the Basque Country (UPV/EHU), 48013 Bilbao, Spain d Basque Center for Applied Mathematics (BCAM), 48009 Bilbao, Spain Keywords: Molecular docking problem Bio-inspired optimization Drug discovery Evolutionary computation Swarm intelligence a b s t r a c t Molecular docking is a Bioinformatics method based on predicting the position and orientation of a small molecule or ligand when it is bound to a target macromolecule. This method can be modeled as an optimization problem where one or more objectives can be defined, typically around an energy scoring function. This paper reviews developments in the field of single- and multi-objective metaheuristics for efficiently addressing molecular docking optimization problems. We comprehensively analyze both problem formulations and applied techniques from Evolutionary Computation and Swarm Intelligence, jointly referred to as Bio-inspired Optimization. Our prospective analysis is supported by an experimental study dealing with a molecular docking problem driven by three conflicting objectives, which is tackled by using different multi-objective heuristics. We conclude that genetic algorithms are the most widely used techniques by far, with a noted increasing prevalence of particle swarm optimization in the last years, being these last techniques particularly adequate when dealing with multi-objective formulations of molecular docking problems. We end this experimental survey by outlining future research paths that should be under target in this vibrant area. 1. Introduction In the last few decades, the field of Bioinformatics has experienced an exponential growth of the amount of biological data available [1], leading to the rise of new challenges related to data storage, management, extraction and analysis. In particular, many Bioinformatics problems (e.g., molecular docking, protein structure prediction, phylogenetic inference, etc.) can be defined in terms of one or more objective functions that have to be minimized or maximized [2]. Most of such problems have been addressed in the literature by resorting to meta-heuristic techniques [3]. In the area of Genomics, the increase of the number of human DNA sequences available has allowed the detection of genes [4], which express proteins that can be considered as therapeutic targets for drug discovery. These candidate proteins are extracted by applying techniques, such as X-crystallography and Nuclear Magnetic Resonance spectroscopy (NMR). This has enabled to know many structural details of ligand–protein complexes, as well as the creation of databases such as the Protein Data Bank (PDB) [5]. This database stores information about biological macromolecule crystal structures, whose atom coordinates are utilized in in silico studies. These advances on laboratory and computational techniques have shed light on the approaches applied in drug discovery, and particularly molecular docking. This procedure consists in finding the best orientation and position of a given chemical compound (ligand) to a therapeutic target (macromolecule). ∗ Correspondence to: P. Tecnologico, Edificio 700, 48160 Derio, Bizkaia, Spain. E-mail address: [email protected] (J. Del Ser). Molecular docking has significantly contributed to structurebased drug discovery as has been emphasized [6], postulated [7] and shown with successful protocols [8] in manifold contributions. Finding the ligand’s best orientation and binding site can be formulated as an optimization problem, whose huge complexity has stimulated the adoption of non-exact search methods belonging to the family of meta-heuristics, such as Evolutionary Algorithms. Indeed, meta-heuristics have been applied to single- and multi-objective formulations of this problem, which is the main focus of this survey. In this regard, addressing the optimization of a molecular docking problem requires to deal with two issues: the definition of the energy scoring function and the choice of the search method. The energy scoring function evaluates the quality of the ligand–macromolecule interactions through approximations to provide a compromise between accuracy and speed. The definition of energy scoring functions has been conducted over the last two decades to yield different main categories [9]: (1) force field-based methods, (2) empirical scoring functions, (3) knowledge-based potentials, and (4) descriptor-based scoring functions. These approaches can model the molecular docking simulations as rigid or flexible macromolecule–ligand conformations only if flexibility can be applied as degrees of freedom in the ligand or macromolecules bonds. Simulations assuming flexibility are more realistic than those in which flexibility is not endowed. The quality of the molecular docking solutions are measured in terms of the energy associated to the ligand– protein interactions (∆G), which indicates the stability of the complex computed. Likewise, the RMSD (root-mean-square deviation) is also employed as a mean to gauge the distance between the co-crystallized ligands and the computed ones, measured in Angstroms (Å). The search method refers to the technique that is used to find the best solution, which represents the best ligand’s conformation to the macromolecule in terms of energy. As anticipated above, the most widely adopted family of solvers for molecular docking are meta-heuristics, which have been used since the 90’s [10]. Nevertheless, there is a lack of studies in which this kind of optimization techniques are reviewed towards providing biologists and researchers in applied metaheuristics with an organized corpus of references related to this research area, including a prospect of different perspectives and challenges. Therefore, our motivation in this paper is to provide the audience with a comprehensive survey of contributed works related to bio-inspired meta-heuristic algorithms applied the molecular docking optimization problem. Our bibliographic study provides an overall overview, which we deem of utmost interest for both biologists and professionals in bioinformatics with strong background in the problem, as well as for computer scientists working on the wide field of bio-inspired meta-heuristics. The first target group of readers will find in this manuscript informed hints on which techniques can be more promising to be applied in practice, whereas the second group will discover a vast spectrum of features related with this real-world optimization problem, that can be useful to investigate and derive new heuristic techniques. The contributions of this paper can be summarized as follows: •The proposed survey focuses on different yet complementary perspectives: the different scoring functions that have been adopted in the literature, the diverse single/multiobjective formulations derived to model this problem, and the variety of bio-inspired algorithms used to deal with them. More than 115 references have been selected for revision, and catalogued accordingly. •We provide practical guidelines to the audience by presenting and discussing results of an experimental case study based on a three-objective formulation of the molecular docking problem on 75 instances with multi-objective metaheuristics. A thorough algorithmic comparison with recommendations is included to support experts in the selection of algorithms. •A prospective discussion is made to elaborate on open issues, trends and challenges with the aim of encouraging research community in this area to follow-up with molecular docking optimization. The remainder of this paper is structured as follows: first, a background section devoted to bio-inspired optimization algorithms is included. Next, different categories of energy scoring functions are introduced, which is used as the baseline criterion to classify the literature in this field until today. Then we introduce and discuss different single- and multi-objective heuristic approaches applied to solve the molecular docking problem by sorting them by the number of objectives optimized in each reported contribution. Once we have examined the state of the art, we present and discuss the results obtained for the aforementioned molecular docking problem with three objectives. Finally, we discuss open issues, trends and challenged, and draw final remarks on the future of this particularly active research area. 2. Bio-inspired optimization algorithms Bio-inspired metaheuristics [11] for solving complex optimization problems have become a very popular research topic in last decades. Among the myriad of methods proposed to date, classical techniques from the wide field of Evolutionary Computation – such as Genetic Algorithms (GA) or Differential Evolution – are arguably the most established solvers because of their versatility and relatively longer presence in the community. However, modern bio-inspired approaches such as Particle Swarm Optimization (PSO) [12], Artificial Bee Colony (ABC) [13], Cuckoo Search (CS) [14] and Firefly Algorithm [15] have started to attract more attention for continuous- and discrete-variable optimization problems, due to their prominent performance shown in manifold real-world problems. Most of these new metaheuristics hinge on search mechanisms that are rather different than those characterizing evolutionary search techniques, as we will expose briefly in what follows. There are three main phases in most bio-inspired evolutionary algorithms: (1) initialization of each individual in a set of candidate solutions to the problem at hand (population), which is randomly generated according to some numerical representation (solution encoding); (2) each solution in the population is evaluated in terms of its fitness value, which can be used for ranking the individual solution within the population for the purpose of selecting or discarding it along the search process, (3) generation of a new population by recombination and/or mutation of solutions in the existing population, followed by the evaluation of the newly produced solutions (offspring) as per the fitness metric of the problem; (4) a replacement criterion is adopted and applied to select which solutions remain in the population and which ones are discarded due to their worse fitness. Selection, recombination, mutation, evaluation and replacement are repeated until a predefined stop condition is met. A general pseudocode of the whole optimization process undertaken by bio-inspired evolutionary algorithms is shown in Algorithm 1. At the end of this process, the best solution found so far along the search process is returned, which is declared to be the solution provided by the algorithm. In regards to newer search heuristics like the ones exemplified above, their search process is rather guided by the interactions Algorithm 1 Pseudocode of a generic evolutionary algorithm. 1: initializePopulation() 2: evaluatePopulation() 3: while not StopCondition() is met do 4: selectParents() 5: recombination() 6: evaluateOffspring() 7: replaceIndividuals() 8: end while 9: returnGlobalBestSolution() between multiple agents endowed with simple behavioral rules for interaction and communication with their counterparts. It is the collective knowledge acquired by a set (swarm) of such agents where the search power of the overall algorithm resides, as their performed exploration benefits from the information about the problem transmitted and conveyed throughout the members of the swarm. A biological simile often lies underneath this wide class of solvers, such as the behavioral patterns observed in bird flocks or fish schools. This gives rise to the concept of Swarm Intelligence [16], which refers to all search algorithms that conform to the multi-agent procedure sketched previously. 3. The molecular docking problem This section is devoted to the molecular docking problem. It is divided into three different subsections, concentrated on: (1) solution encoding strategy; (2) the energy scoring function, which can be classified according to four main categories; (3) the objectives to be minimized; and (4) the search method used to find the best ligand–macromolecule interactions (solutions). 3.1. Solution encoding Molecular docking processes are aimed at finding 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 side chains of the receptor. Fig. 1 illustrates the most commonly used strategy for solution encoding, which consists of a real-value vector of 7+nvariables. In this vector, 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. In order to reduce the computational cost, a grid-based methodology can be carried out by which the protein active site is embedded in a 3D rectangular grid, and by taking into account all the protein atoms. In this way, the protein contribution at any given point is obtained by trilinear interpolation in the grid cell. As a consequence of this process, the range of translation variables (x,y,z) is [0· · · 120], which bounds effectively the limits of the coordinates of the grid space previously set for each problem (and in essence, reduces the size of its search space). Orientation (quaternion) and torsion variables are measured in radians and encoded in the range of [−π, π]. 3.2. Energy scoring functions for the molecular docking problem The energy scoring function allows evaluating the solutions to the problem (ligand–macromolecule interactions), making approximations that establish a trade-off between accuracy and speed. In this context, a complete review by Liu et al. [9] categorizes different energy scoring functions reflecting the current progression in this field. This categorization is composed by the following categories: (1) force-field-based energy functions, (2) empirical or regression-based energy functions, (3) knowledgebased scoring functions and (4) descriptor-based scoring functions. Functions that were first developed are included in the category of force-field-based energy functions proposed by Martin Karplus et al. [17,18]. This type of energy scoring function was proposed to calculate the ligand–macromolecule interactions through the inclusion of non-covalent terms, like the so-called Van der Waals interactions, the electrostatic energy and H-bond terms. Some molecular docking software programs incorporating this evaluation function are AutoDock [19], DOCK [20], COMBINE [21], GoldScore [22], MedusaScore [23] and CHARMM [24]. The following general functional form became more widely adopted for this type of energy function: ∆G=∆Evdw+∆Eelectrostatic +∆EH−bond +∆Edesolvatation (1) The second type of evaluation energy functions is referred to as empirical scoring function. There are a few of exemplifying software frameworks equipped with this class of scoring functions, such as ChemScore [25], X-Score [26] and GlideScore [27]. The total score is composed of a sum of scores based on H- bond interactions, bonds in which metals and lipophilic atom contacts are involved and a set of penalty scores based on rigid rotatable bonds, ligand–macromolecule steric clashes, and other terms based on ligand–macromolecule interactions’ constraints that could be useful, for example, to limit the search space of the considered search method (see Section 3.4). An example of equation that represents this method is as follows: Stotal =SH−bond +Smetal +Slipophilic +Protatable +Pclashes (2) A third category of energy functions corresponds to knowledge-based scoring function, as they are known in the reviewed literature. Software packages using these functions include DrugScore [28] and IT-SCORE [29]. The value of this energy function builds upon the sum of pairwise statistical potentials between protein and ligand, e.g.: A= lig ∑ n macro ∑ m ωij(r),(3) where lig and macro refer to the ligand and the macromolecule, respectively. A fourth category is based on machine learning based-scoring functions that use modern Quantitative Structure–Activity Relationship (QSAR) analysis into the ligand–protein interactions. This type of energy function extracts and exploits descriptors from the ligand and macromolecule to improve the solutions of molecular docking simulations. This is an emerging strand in the area that can be very promising if combined with other energy scoring functions and/or search methods, such as meta-heuristics. 3.3. Objectives of the molecular docking problem The main objective of the molecular docking problem is to find the best conformation between a ligand (L) and a receptor (R) in terms of energy and RMSD. This goal can be formulated as a single- or a multi-objective optimization problem. Therefore, the 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. molecular docking problem can be formally expressed by letting A and Bdenote the ligand and protein molecule, respectively. If we define f:R→Ras a scoring function that evaluates solutions with regards to a given score, and let Cbe the set of all feasible molecular docking solutions between Aand B, then the molecular docking problem aims at finding a conformation −→ x∈Cthat satisfies f(−→ x)≤f(−→ y)∀−→ y∈C:−→ y= −→ x. According to the reviewed literature, if molecular docking is formulated as a single-objective optimization problem, the fitness function to optimize is usually set to the final binding energy ∆G. However, in the middle of 2000’s, Grosdidier et al. [30] and Janson et al. [31] proposed new approaches in which multiobjective optimization was applied to the problem at hand by simultaneously minimizing several objectives that will be described below. These studies were subsequently followed by the work performed by Liu et al. [32] and Boisson et al. [33]. There are other remarkable studies published in the last decade that embraced this design choice, such as [12,34–37]. More recently, a new approach to the problem was proposed by López-Camacho et al. [38], where three objectives were first considered within the formulated molecular docking problem, thereby unleashing new perspectives over the field of multi-objective optimization strategies that could potentially improve the accuracy of drug discovery. In a general multi-objective problem, it is necessary to define two spaces: (1) the decision space, that includes the set of feasible solutions; and (2) the objective space, which involves the fitness (scoring) values of the objectives to be optimized. Objectives formulated for molecular docking in the literature can be enumerated as follows: •In Grosidier et al. [30], objectives that were minimized correspond to two scores resulting from a SimpleFitness and aFullFitness. The first objective computes the sum of the intramolecular energy (Eintra) of Rand L, and the intermolecular energy (Einter ) of the complex R-Lin Expression (4). The second objective is a term that corresponds to the solvation energy ∆Ggiven in Expression (5): SimpleFitness =ER intra +EL intra +Einter ,(4) FullFitness =SimpleFitness +∆Gelec.solv.(5) •In Janson et al. [31], Sandoval-Perez et al. [34] and García- Godoy et al. [36], the objectives used were: Eintra =EL intra+R intra,(6) Einter =ER−L intra+R−L intra,(7) where Eintra refers to the unbound and bound states of the R and L, and Einter refers to the unbound and bound states of the complex R-L. •In Boisson et al. [33], the authors proposed to optimize Einter as per (7), and a surface criterion based on the hypothesis that the surface of Land Rdecreases because of the penetration of Linto R. •Liu et al. [32] elaborated on an approach that optimized the energy associated to the Van der Waals and torsion bonds, namely: EL−R=∆Sconf +WvdW ,(8) and a term based on the geometry dissimilarity (GD): GS =√∑N i=1∥xindiv,i−xinput,i∥2 N,(9) where xindiv,iand xinput,iare the positions of the ith heavy atom in each individual and the input L-Rconformer. •Gu et al. [35]: The authors proposed an approach that minimizes the resulting score from three different scoring energy functions (see Section 3.2): the force-field-based, empiricalbased and knowledge-based functions (see Eqs. (1),(2),(3)), which are treated as objectives during the docking optimization. •In López-Camacho et al. [37] and García-Nieto et al. [12] the molecular docking problem is defined on the basis of two different metrics to be optimized, Einter and the RMSD score. This score measures the average distance between the known ligand position in the receptor and the computed position of the docking ligand, which takes into account symmetry, partial symmetry (e.g. symmetry within a rotatable branch) and near-symmetry in a simple heuristic fashion. This score is given by: RMSDab =max(RMSD′ ab,RMSD′ ba),(10) where RMSD′ ab =√1 N∑iminjrij 2, and 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. •In García-Godoy et al. [38], a three-objective molecular docking problem is first posed and solved, based on the simultaneous optimization of Einter ,Eintra and the RMSD score defined above. Interestingly, this preliminary work unveiled the inherent benefits springing from considering a third objective in the molecular docking problem, thus settling the intuition in which this manuscript finds its motivation. 3.4. Meta-heuristics for the molecular docking problem We have reviewed the most relevant works attending the applied meta-heuristics by considering two main categories: single- and multi-objective algorithms. This categorization can be very useful not only for practitioners that are interested in using one approach or another, but also for developers willing to implement innovative search methods. Furthermore, these approaches have been also organized taking into account the standard categorization of scoring functions provided by Liu et al. [9]. 3.4.1. Single-objective algorithms The use of molecular docking tools to predict the ligandinteractions were pioneered in the early 80’s, but it was not until the 90’s when researchers started to apply meta-heuristics to optimize single-objective formulations of this problem, mainly considering the free binding energy ∆G. During this decade, most published works were based on assessing new search methods that used different scoring functions, being one of the most widely used the one defined by Martin Karplus et al. [17,18], who proposed the use of force fields to compute the interactions between ligands and macromolecules. Some of these functions accommodate some degrees of freedom in the ligands’ and macromolecules’ bonds, which increases the complexity of the problem and the computation time required to solve it to a given accuracy, but yields more realistic molecular docking simulations. In this context, studies such as [19,22,39–44] gravitated on the application of GA by minimizing the final binding energy. Previously, we have also mentioned that some of these early force-field scoring functions allow to apply flexibility to the ligands’ conformations as the case of the studies performed in [19,22,39,42,43]. For example, the authors in [39] were the first to apply flexibility to three instances of the molecular docking problem, and the studies provided by Morris et al. [19] questioned which search method is more efficient in escaping from a local minimum, concluding that the Lamarckian Genetic Algorithm (LGA), a GA with Local Search (LS) obtained the best solutions, followed by GA and Simulated Annealing (SA), which failed to return better solutions in terms of energy and RMSDs. This decade is also characterized by studies based on the use of the Monte Carlo algorithm. For example, [45] used a forcefield based fitness function and Monte Carlo as search method in flexible ligands’ instances. In [46], the authors also proposed a Monte Carlo algorithm to predict 20 molecules’ interactions where flexibility was assumed. In [47], a comparative analysis was carried out to assess the performance of several search methods: GA, Tabu Search (TS), SA and an Evolutionary Programming (EP) algorithm. The authors showed that the GA outperforms the rest of algorithms in terms of energy. Furthermore, the provided discussion suggested that GA is a very effective local search algorithm, but prone to falling in low-energy local minimum. TS, on its own, seems to perform a more global search, which permits to locate a global energy minimum, but less accurately than GA. The authors noted the importance of publishing studies about the search methods to improve future techniques to be applied to molecular docking simulations. In the early 2000’s, the use of EAs continued playing a prominent role, and a tendency of performing comparative analysis emerged at that time. For example, in [48], the authors proposed an EA that used global and local search strategies by integrating decreasing mutations and self-adaptative mutations. This approach was applied to a use case based on a flexible ligand-rigid protein complex. In [49], an EA was designed to minimize Einter , evaluating the solutions by means of the CHARMM tool. The method is applied to three instances of the problem. Budin et al. [50] applied a method based on an EA by optimizing the energy of H-bonds of a set of complexes where flexibility was applied to the ligand (10 rotatable bonds). In [51], a new optimizer was proposed based on GA, Global Underestimator Method and Pole and applied to a set of 27 rigid instances. Solutions were evaluated by a force-field-based energy scoring function. In [28], the authors revolved on a new scoring function, which is part of the knowledge-based scoring function group, as well as on the LGA solver provided by AutoDock. Other remarkable studies include [52,53]. The former uses a multi-population GA and a force-field energy scoring function to evaluate the solutions from a set of 11 complexes with flexible ligands. The latter applies an GA-like method with an empirical energy function for seven complexes, where flexibility is applied to the ligand. At this point, we can conclude that these published works have several features in common: (1) the number of complexes is not enough to do a comparative analysis, (2) flexibility is only applied to the ligand’s conformation; and (3) the number of algorithms applied to perform the comparative analysis is limited. More recently, there is an uprising trend in the field to use Swarm Intelligence techniques rather than Evolutionary Computation approaches to tackle molecular docking based on optimizing a single objective. For example, Nesamalar and Chandran [13] introduced a Genetic Clustering combined with Bee Colony Optimization (BCO). The approach, coined as (GCBCO), was tested by using 10 docking instances from the PDB bind core set. The results showed that GCBCO outperforms other Swarm Intelligence algorithms (specifically, PSO and Ant Colony Optimization - ACO) in terms of energy. The authors in [54] utilized the AutoDock 4.2 energy function to include flexibility in the macromolecule and ligand. Another novel ingredient of this recent reference was a new search algorithm hybridizing features of a quantumbehaved particle swarm optimization (QPSO) and the local search method of Solis and Wets. A thorough analysis of the algorithm’s performance was made by comparing the results obtained by QPSO with those rendered by the LGA optimizer embedded in AutoDock. 3.4.2. Multi-objective algorithms The first studies to consider a multi-objective optimization approach appeared in 2007 [30,31]. In [30], Grosididier et al. presented EADock, which is a multi-objective formulation based on an EA and the minimization of two objectives from two different fitness functions: the total binding energy and the solvation energy as per Eqs. (4) and (5), respectively. In [31], a multiobjective approach based on PSO (namely, ClutMPSO) was shown to minimize efficiently Einter and Eintra. The ligand–macromolecule interactions were evaluated with the early energy function of AutoDock 3.0. Other notable contributions to the multi-objective formulation of the molecular docking problem are [34,36], where the same objectives as in [31] are optimized, and [35]. In the first work, a force-field-based energy function is used to evaluate solutions, and the NSGA-II multi-objective meta-heuristic is used as the search method. In the second work, 11 instances were tested by multi-objective solvers assuming flexibility in both the ligand and macromolecule. Search methods that were used ranged from NSGA-II to ssNSGA-II, SMPSO, MOEA/D and SMS-EMOA. Interestingly SMPSO, a Swarm Intelligence multi-objective optimization method that resembles the search mechanisms underlying the classical PSO solver, was shown to attain outperforming results in terms of energy and RMSD, when compared with the rest of algorithms. In the third work, the authors proposed a multiobjective approach (MoDOCK) that minimizes three fitness scores from different scoring functions. This approach is compared with single-objective approaches that used popular scoring functions, such as: DOCK, Glide, SurFlex, DOCK6 and DOCK6-F. The analysis indicated that MoDOCK dominated the rest of approaches in terms of RMSD. In [37,55], the authors propose a multi-objective approach to optimize Einter and the RMSD score using the energy function of AutoDock 4.2 to enable flexibility in the macromolecule and the ligand. The selected algorithm were NSGA-II, SMPSO, GDE3, and MOEA/D. In [37], it was again observed that SMPSO showed the best overall results in terms of energy and RMSD (value lower than 2 for successful docking results). Finally, in the most recently study reported in [38] three-objectives were considered for minimization: Einter ,Eintra and RMSD. The solvers utilized in this late study were SMPSO, MOEA/D and MPSO/D, whereas the AutoDock 4.2 function was considered to model flexibility in the macromolecule and the ligands. The obtained results in this preliminary study laid the foundations for more elaborated experiments as the one later presented and discussed in this manuscript. 4. A survey of molecular docking using bio-inspired optimization This section classifies the most relevant articles that deal with the molecular docking problem using bio-inspired metaheuristics, by conforming to the topics developed in the previous sections. In this way, the literature can be classified depending on three different taxonomic aspects: the first one deals with the scoring function used to evaluate the solutions produced by the different bio-inspired meta-heuristics; the second one reduces to the particular problem objectives that are under consideration in each work. The third classification criterion takes into account which family of bio-inspired meta-heuristic techniques have been applied to solve molecular docking problem. The main source of information we have used to find the papers considered in this survey has been ScopusR .1The search terms utilized to query this database have been composed by combinations of ‘‘molecular docking’’, ‘‘protein docking’’, ‘‘genetic algorithm’’, ‘‘evolutionary algorithm’’, ‘‘swarm intelligence’’. From all the retrieved results, we have selected the most representative works for reaching informed insights and conclusions about the general trends of this literature strand. As a starting point of our bibliographic analysis, Table 1 organizes related work according to scoring functions that have been used along the years. These scoring functions can be divided into four types as discussed in Section 3.2: force-field-based energy functions, empirical or regression-based energy functions, knowledge-based scoring functions and descriptor-based scoring functions. Works using force-field based energy functions are more common in the area, being the AutoDock scoring function the most prominent among them. Some studies resort to their own formulation of the force-field based function. The bibliographic taxonomy follows in Table 2, although in this case the literature comprising multi-objective approaches 1Scopus: https://www.scopus.com. is categorized in terms of the selected optimization objectives, bearing in mind that when one single objective is considered, the intermolecular energy Einter is always chosen. Table 3 inspects the use of different bio-inspired metaheuristics to solve the molecular docking problem from the 90’s to the present day. It is remarkable that the use of GA has been prominent all over the history of the field, arguably by the particularly good results achieved by this solver when addressing molecular docking problems. Some approaches using the Monte Carlo algorithm (MC) were noticed at the dawn of the field, but this approach has fallen into disuse. PSO techniques have been used more frequently in recent years, as they were empirically found to improve existing methods at the time. It is worth noting that LS approaches have been also hybridized with other common meta-heuristics, hence tracing a promising research path around Memetic Computing [104], a paradigm within Bio-inspired Optimization that precisely embraces this hybridization as its core design principle. 5. Experimental study: optimizing three objectives In the previous section, we have described the historical evolution of the research activity around the use of bio-inspired meta-heuristics to deal with the molecular docking problem. To provide readers interested in this subject with empirical evidences of the practical benefits derived from this research area, in this section we conduct a study where three state-of-the-art bioinspired meta-heuristics are used to solve a molecular docking problem considering three objectives: the inter-molecular energy Einter , the intra-molecular energy Eintra and the Root Median Square Deviation (RMSD) score. As commented previously, Einter relates to the ligand–receptor affinity, Eintra characterizes the ligand deformity, and the RMSD score measures the difference of atomic distances between the co-crystallized ligand and the computed one. The study involves both, a quantitative analysis of the solutions found by the algorithms (computer scientist perspective) and a qualitative analysis of some selected results (biologist perspective). Three multi-objective algorithms were selected for the experimental part of this study due to their renowned good performance in other multi-objective problems: SMPSO and MOEA/D, which have yielded good results when dealing with this problem in the recent past [36]; and MPSO/D, recently applied to the same problem optimizing only two objectives (Einter and RMSD) and using rigid ligand–protein complexes [12]. We next describe them briefly for the sake of completeness: •SMPSO is a multi-objective PSO algorithm that incorporates a velocity restriction mechanism. The approach allows producing new effective particle positions in those cases where the velocity becomes too high. This solver also includes additional features, such as the use of polynomial mutation as a turbulence factor and an external archive to store non-dominated solutions found during the search. In [116], a comparative analysis of SMPSO with five representative multi-objective solvers was performed in terms of the quality of the resulting approximation sets and the convergence speed to the Pareto front. The obtained results showed that SMPSO returned the best overall results in such terms. •MOEA/D is a multi-objective EA based on decomposition [117]. This algorithm divides the multi-objective problem into a number of scalar optimization subproblems, which are solved simultaneously by evolving the population. At each generation, the best solution for each subproblem is obtained. The neighborhood relations among these subproblems are defined based on the distances between their aggregation coefficient vectors. The optimization of each subproblem is based on the information from the subproblems within its neighborhood. Table 1 Classification of the literature in terms of scoring functions. Year Reference Force field Empirical Knowledge (DrugScore) Descriptor (QSAR) Other DOCK GOLD AutoDock AMBER GROMOS 1994 Xiao and Williams [56] Xiao and Williams [57] Knegtel1 et al. [58] 1995 Oshiro et al. [39] Clark and Ajay [59] Jones et al. [60] Gehlhaar et al. [42] Read et al. [46] Verkhivker et al. [41] 1997 Westhead et al. [47] Maddalena and Snowdon [61] Jones et al. [22] Levine et al. [44] 1998 Vieth et al. [62] Morris et al. [19] 1999 Chen and Chi [63] Liu and Wang [45] 2000 Yang and Kao [48] Taylor and Burnett [49] 2001 Budin et al. [50] David et al. [51] Pegg et al. [64] 2002 Hou et al. [65] Sotriffer et al. [28] 2003 McGann et al. [66] 2004 Lee et al. [67] Li et al. [52] Yang and Chen [53] Magalhães et al. [68] 2005 Li et al. [69] Janson and Merkle [70] Li et al. [71] Li et al. [72] Chang et al. [73] 2006 Wiley et al. [74] Oduguwa et al. [75] Chen et al. [76] 2007 Grosdidier et al. [30] Zhao and Sanner [77] Kroemer [78] 2008 Kang et al. [79] Sung [80] Janson et al. [31] 2009 Liu et al. [32] Thiriot and Monard [81] Liu et al. [82] Tavares et al. [83] 2010 Fuhrmann et al. [84] Ling et al. [85] Atilgan and Hu [86] Archetti et al. [87] Meier et al. [88] Boisson et al. [33] Mesmoudi et al. [89] 2011 Rondón et al. [90] 2012 Liu et al. [91] Lima et al. [92] Nesamalar and Chandran [13] Koohi-Moghadam and Rahmani [93] Liu et al. [94] 2014 Zhen-yu et al. [95] 2015 Fu et al. [54] Li et al. [96] Ravindranath et al. [97] Xu et al. [98] García-Godoy et al. [36] López-Camacho et al. [99] 2016 López-Camacho et al. [37] García-Nieto et al. [55] Guan et al. [100] 2018 Leonhart et al. [101] Guan et al. [102] García-Nieto et al. [12] Lopez-Camacho et al. [38] Fu et al. [103] Table 2 Classification of the literature in terms of optimization objectives. Year Reference Eintra Einter ∆elec solvSL−R∆Sconf EvdW GS RMSD 2007 Grosidier et al. [30]  2008 Janson et al. [31]  2010 Boisson et al. [33]  2013 Sandoval-Perez et al. [34]  2015 García-Godoy et al. [36]  Liu et al. [9]   Gu et al. [35] N/A N/A N/A N/A N/A N/A N/A N/A 2016 López-Camacho et al. [37]  2018 García-Nieto et al. [12]  Table 3 Classification of the literature in terms of the utilized meta-heuristics (ABC: Artificial Bee Colony) Year Reference MC EA TS SA PSO ACO ABC LS GA DE Other 1994 Xiao and Williams [56] Xiao and Williams [57] Knegtel1 et al. [58] 1995 Oshiro et al. [39] Clark and Ajay [59] Willett [40] Jones et al. [60] Gehlhaar et al. [42] Read et al. [46] 1996 Clark and Westhead [43] Verkhivker et al. [41] 1997 Westhead et al. [47]   Maddalena and Snowdon [61] Jones et al. [22] Levine et al. [44] 1998 Vieth et al. [62]  Morris et al. [19]   Szczerbicka et al. [105] 1999 Chen and Chi [63]  Liu and Wang [45] 2000 Yang and Kao [48] Taylor and Burnett [49] 2001 Budin et al. [50] David et al. [51]  Pegg et al. [64] 2002 Hou et al. [65]   Sotriffer et al. [28]  2004 Lee et al. [67]  Schneidman-Duhovny et al. [106] Magalhães et al. [68] 2005 Li et al. [69] Janson and Merkle [70]   Li et al. [71] Li et al. [72] Chang et al. [73]  2006 Wiley et al. [74]  Oduguwa et al. [75] Chen et al. [76]  2007 Grosdidier et al. [30] Zhao and Sanner [77] Kroemer [78]   2008 Kang et al. [79] Sung [80] Janson et al. [31] 2009 Liu et al. [32] Thiriot and Monard [81] Liu et al. [82] Tavares et al. [83]   2010 Fuhrmann et al. [84] Ling et al. [85] Atilgan and Hu [86] Archetti et al. [87] Chang et al. [107] Meier et al. [88] Boisson et al. [33] Mesmoudi et al. [89] (continued on next page) Table 3 (continued). Year Reference MC EA TS SA PSO ACO ABC LS GA DE Other 2011 Rondón et al. [90] 2012 Lima et al. [92] Nesamalar and Chandran [13]   Liu et al. [94]   2013 Sandoval-Perez et al. [34] López-Camacho et al. [108]   Rakshit et al. [109]   Liu et al. [110] 2014 Magalhaes et al. [111] Zhen-yu et al. [95] 2015 Fu et al. [54]  Li et al. [96] Ravindranath et al. [97] Xu et al. [98] García-Godoy et al. [36]    López-Camacho et al. [99]   2016 López-Camacho et al. [37]   García-Nieto et al. [55] Guan et al. [100] García-Godoy et al. [112]   Peh and Hong [113] 2017 Cecilia et al. [114] Guan et al. [115] 2018 Leonhart et al. [101] Guan et al. [102]  García-Nieto et al. [12] Lopez-Camacho et al. [38]  Fu et al. [103] •MPSO/D is a relatively new multi-objective PSO algorithm based on decomposition [118]. This algorithm decomposes the objective problem into a set of sub-regions on the basis of a set of direction vectors. Each subregion has a solution that maintains a level of diversity. MPSO/D also exploits the crowding distance to calculate the fitness of the solutions for the selection operator, as well as the neighboring particles (solutions) of a given particle to infer the global best historical position found by the algorithm. These three multi-objective algorithms have been implemented in C++ as part of the jMetalCpp2project [108], which is essentially a software framework for single- and multiobjective optimization meta-heuristics. The energy score evaluation has been performed by AutoDock 4.2 [119] which, as mentioned before, is a widely used tool for virtual drug discovery involving rigid and flexible docking simulations. The integration of AutoDock with jMetalCpp3provides a flexible capacity to use the majority of the meta-heuristics implemented for the molecular docking problem. Specifically, jMetalCpp is in charge of evolving solutions as per the search operators of the metaheuristic algorithm under choice, only delegating in AutoDock the evaluation of a newly produced candidate. This experimental use case was carried out over an extensive set of 75 molecular complexes that have been used in many other studies [119]. Table 4 summarizes this set of instances, including the PDB codes of these structures and the range of the crystalographic resolution in Å. Moreover, these complexes are classified in terms of the type of their ligand: small,medium, and large size inhibitors, as well as cyclic urea inhibitors. All these compounds are characterized by a flexibility of the ligand limited to 10 torsional degrees of freedom. These instances were previously preprocessed (both the macromolecule and the ligand), and their potentials were calculated for a three-dimensional square 2http://jmetalcpp.sourceforge.net/ (accessed on November 15th, 2018). 3http://khaos.uma.es/autodockjmetal/ (accessed on November 15th, 2018). (x,y,z) grid of 120 Å per dimension and a grid spacing of 0.375 Å, using the AutoGrid4 software. Once all the macromolecules and ligands have been prepared, we executed 30 independent runs for each pair of algorithm and instance. Table 5 indicates the parameter values of the three algorithms included in the benchmark. Each used a population/swarm size of 100 individuals (or particles), and were limited to a maximum of 1,000,000 fitness evaluations. The rest of parameters were set with their nominal values. We used two quality indicators to evaluate the algorithms’ performance: the Hypervolume (IHV ) [120] and the Unary Additive Epsilon Indicator (Iϵ+) [121]. IHV gives a measure of both convergence and diversity, whereas the Iϵ+takes into account the convergence degree of the obtained Pareto front approximations. As molecular docking is a real-world problem, the Pareto front is not known a priori. Consequently, in order to calculate these two measures for each instance, a reference front was generated using all the non-dominated solutions obtained for an instance using all executions of the three algorithms. For these indicators, we calculated the median and the Inter-Quartile Range (IQR) as statistical measures of central tendency and dispersion. Tables 6 and 7list such statistics for IHV and Iϵ+, respectively. The cells shaded in dark gray indicate the algorithm scoring the best median value for the corresponding indicator, whereas those in light gray refer to the algorithm achieving the second best median indicator value for each instance. In the IHV indicator, the higher the convergence and diversity degree of a front is, the higher (better) the resulting IHV value will be. We can observe in Table 6 that MOEA/D obtained the best median value in 17 instances and the second best value in 18 instances. SMPSO obtained the best and second best median value in 38 and 17 instances, respectively. Finally, MPSO/D achieved the best and second best median value 20 and 40 times over the benchmark. Taking into account the Iϵ+indicator, the higher the convergence of a front is, the lower (and better) the value of Iϵ+will be. Having this in mind, it can be observed that MOEA/D obtained 14 best median values and 18 s best median values of this indicator. SMPSO obtained the better results with 43 and [88] R. Meier, M. Pippel, F. Brandt, W. Sippl, C. Baldauf, ParaDockS: a framework for molecular docking with population-based metaheuristics, J. Chem. Inform. Model. 50 (5) (2010) 879–889, http://dx.doi.org/10.1021/ ci900467x. [89] S. Mesmoudi, J. Tavares, L. Jourdan, E.-G. Talbi, Variable genetic operator search for the molecular docking problem, in: Evolutionary Computation, Machine Learning and Data Mining in Bioinformatics, Springer Berlin Heidelberg, 2010, pp. 1–12, http://dx.doi.org/10.1007/978-3-642-12211- 8_1. [90] P. Rondon, H. Arguello, R. Torres, Development of a zoned genetic algorithm for semi-flexible protein-ligand docking in drug design, in: 2011 6th Colombian Computing Congress, CCC, IEEE, 2011, pp. 1–4, http://dx.doi.org/10.1109/colomcc.2011.5936321. [91] Y. Liu, L. Zhao, W. Li, D. Zhao, M. Song, Y. Yang, FIPSDock: A New molecular docking technique driven by fully informed swarm optimization algorithm, J. Comput. Chem. 34 (1) (2012) 67–75, http://dx.doi.org/10. 1002/jcc.23108. [92] A.N. Lima, E.A. Philot, D. Perahia, A.S.K. Braz, L.P. Scott, GANM: A Protein– ligand docking approach based on genetic algorithm and normal modes, Appl. Math. Comput. 219 (2) (2012) 511–520, http://dx.doi.org/10.1016/ j.amc.2012.06.030. [93] M. Koohi-Moghadam, A.T. Rahmani, Molecular docking with oppositionbased differential evolution, in: Proceedings of the 27th Annual ACM Symposium on Applied Computing - SAC’12, ACM Press, 2012, pp. 1387–1392, http://dx.doi.org/10.1145/2245276.2231996. [94] Y. Liu, W. Li, R. Ma, Particle swarm optimization on flexible docking, Int. J. Biomath. 05 (05) (2012) 1250044, http://dx.doi.org/10.1142/ s1793524511001866. [95] L. Zhen-yu, C. Ming-hua, K. Ling, Application of immune algorithm in molecular docking, in: Proceedings 2014 IEEE International Conference on Security, Pattern Analysis, and Cybernetics, SPAC, IEEE, 2014, pp. 122–125, http://dx.doi.org/10.1109/spac.2014.6982670. [96] Z. Li, J. Gu, H. Zhuang, L. Kang, X. Zhao, Q. Guo, Adaptive molecular docking method based on information entropy genetic algorithm, Appl. Soft Comput. 26 (2015) 299–302, http://dx.doi.org/10.1016/j.asoc.2014.10. 008. [97] P.A. Ravindranath, S. Forli, D.S. Goodsell, A.J. Olson, M.F. Sanner, AutoDockFR: advances in protein-ligand docking with explicitly specified binding site flexibility, in: J.S. Fetrow (Ed.), PLOS Comput. Biol. 11 (12) (2015) e1004586, http://dx.doi.org/10.1371/journal.pcbi.1004586. [98] W. Xu, A.J. Lucke, D.P. Fairlie, Comparing sixteen scoring functions for predicting biological activities of ligands for protein targets, J. Mol. Graph. Modelling 57 (2015) 76–88, http://dx.doi.org/10.1016/j.jmgm.2015.01. 009. [99] E. López-Camacho, M.J.G. Godoy, J. García-Nieto, A.J. Nebro, J.F. Aldana- Montes, Solving molecular flexible docking problems with metaheuristics: A comparative study, Appl. Soft Comput. 28 (2015) 379–393, http://dx. doi.org/10.1016/j.asoc.2014.10.049. [100] B. Guan, C. Zhang, J. Ning, EDGA: a population evolution direction-guided genetic algorithm for protein–ligand docking, J. Comput. Biol. 23 (7) (2016) 585–596, http://dx.doi.org/10.1089/cmb.2015.0190. [101] P.F. Leonhart, E. Spieler, R. Ligabue-Braun, M. Dorn, A biased random key genetic algorithm for the protein–ligand docking problem, Soft Comput. (2018) http://dx.doi.org/10.1007/s00500-018-3065-5. [102] B. Guan, C. Zhang, Y. Zhao, An efficient ABC_DE_based hybrid algorithm for protein–ligand docking, Int. J. Mol. Sci. 19 (4) (2018) 1181, http: //dx.doi.org/10.3390/ijms19041181. [103] Y. Fu, Z. Chen, J. Sun, Random drift particle swarm optimisation algorithm for highly flexible protein-ligand docking, J. Theoret. Biol. 457 (2018) 180–189, http://dx.doi.org/10.1016/j.jtbi.2018.08.034. [104] F. Neri, C. Cotta, Memetic algorithms and memetic computing optimization: A literature review, Swarm Evol. Comput. 2 (2012) 1–14, http: //dx.doi.org/10.1016/j.swevo.2011.11.003. [105] H. Szczerbicka, M. Becker, M. Syrjakow, Genetic algorithms: a tool for modelling, simulation, and optimization of complex systems, Cybern. Syst. 29 (7) (1998) 639–659, http://dx.doi.org/10.1080/019697298125461. [106] B.S.P.D. Schneidman-Duhovny, B.S.P.R. Nussinov, B.S.P.H.J. Wolfson, Predicting molecular interactions in silico: II. Protein-protein and proteindrug docking, Current Med. Chem. 11 (1) (2004) 91–107, http://dx.doi. org/10.2174/0929867043456223. [107] D.T.-H. Chang, J.-H. Lin, C.-H. Hsieh, Y.-J. Oyang, On the design of optimization algorithms for prediction of molecular interactions, Int. J. Artif. Intell. Tools 19 (03) (2010) 267–280, http://dx.doi.org/10.1142/ s0218213010000182. [108] E. López-Camacho, M.J. García-Godoy, A.J. Nebro, J.F. Aldana-Montes, jMetalCpp: optimizing molecular docking problems with a C++ Metaheuristic framework, Bioinformatics 30 (3) (2013) 437–438, http://dx.doi. org/10.1093/bioinformatics/btt679. [109] P. Rakshit, A. Konar, A. Chowdhury, E. Kim, A.K. Nagar, Muti-objective evolutionary approach of ligand design for protein-ligand docking problem, in: 2013 IEEE Congress on Evolutionary Computation, IEEE, 2013, pp. 237–244, http://dx.doi.org/10.1109/cec.2013.6557576. [110] Y. Liu, L. Zhao, M. Li, C. Zhao, Swarm intelligence for molecular docking, Int. J. Model. Ident. Control 18 (4) (2013) 357, http://dx.doi.org/10.1504/ ijmic.2013.053541. [111] C.S. de Magalhães, D.M. Almeida, H.J.C. Barbosa, L.E. Dardenne, A dynamic niching genetic algorithm strategy for docking highly flexible ligands, Inform. Sci. 289 (2014) 206–224, http://dx.doi.org/10.1016/j.ins.2014.08. 002. [112] M. García-Godoy, E. López-Camacho, J. García-Nieto, A. Nebro, J. Aldana- Montes, Molecular docking optimization in the context of multi-drug resistant and sensitive EGFR mutants, Molecules 21 (11) (2016) 1575, http://dx.doi.org/10.3390/molecules21111575. [113] S.C.W. Peh, J.L. Hong, Protein ligand docking using simulated jumping, in: Computational Science and Its Applications – ICCSA 2016, Springer International Publishing, 2016, pp. 1–10, http://dx.doi.org/10.1007/978- 3-319-42111-7_1. [114] J.M. Cecilia, J.-M. Cutillas-Lozano, D. Giménez, B. Imbernón, Exploiting multilevel parallelism on a many-core system for the application of hyperheuristics to a molecular docking problem, J. Supercomput. 74 (5) (2017) 1803–1814, http://dx.doi.org/10.1007/s11227-017-1989-7. [115] B. Guan, C. Zhang, J. Ning, Genetic algorithm with a crossover elitist preservation mechanism for protein–ligand docking, AMB Express 7 (1) (2017) http://dx.doi.org/10.1186/s13568-017-0476-0. [116] A. Nebro, J. Durillo, J. Garcia-Nieto, C.C. Coello, F. Luna, E. Alba, SMPSO: A new PSO-Based metaheuristic for multi-objective optimization, in: 2009 IEEE Symposium on Computational Intelligence in Milti-Criteria Decision-Making, IEEE, 2009, pp. 66–73, http://dx.doi.org/10.1109/mcdm. 2009.4938830. [117] Y. Zhao, H.-L. Liu, Multi-objective particle swarm optimization algorithm based on population decomposition, in: Intelligent Data Engineering and Automated Learning – IDEAL 2013, Springer Berlin Heidelberg, 2013, pp. 463–470, http://dx.doi.org/10.1007/978-3-642-41278-3_56. [118] C. Dai, Y. Wang, M. Ye, A new multi-objective particle swarm optimization algorithm based on decomposition, Inform. Sci. 325 (2015) 541–557, http://dx.doi.org/10.1016/j.ins.2015.07.018. [119] G.M. Morris, R. Huey, W. Lindstrom, M.F. Sanner, R.K. Belew, D.S. Goodsell, A.J. Olson, AutoDock4 and AutoDockTools4: Automated Docking with selective receptor flexibility, J. Comput. Chem. 30 (16) (2009) 2785–2791, http://dx.doi.org/10.1002/jcc.21256. [120] E. Zitzler, L. Thiele, Multiobjective evolutionary algorithms: a comparative case study and the strength Pareto approach, IEEE Trans. Evol. Comput. 3 (4) (1999) 257–271, http://dx.doi.org/10.1109/4235.797969. [121] E. Zitzler, L. Thiele, M. Laumanns, C. Fonseca, V. da Fonseca, Performance assessment of multiobjective optimizers: an analysis and review, IEEE Trans. Evol. Comput. 7 (2) (2003) 117–132, http://dx.doi.org/10.1109/tevc. 2003.810758. [122] G. Ramachandran, C. Ramakrishnan, V. Sasisekharan, Stereochemistry of polypeptide chain configurations, J. Mol. Biol. 7 (1) (1963) 95–99, http: //dx.doi.org/10.1016/S0022-2836(63)80023-6,http://www.sciencedirect. com/science/article/pii/S0022283663800236. [123] A. Eiben, S. Smit, Parameter tuning for configuring and analyzing evolutionary algorithms, Swarm Evol. Comput. 1 (1) (2011) 19–31, http: //dx.doi.org/10.1016/j.swevo.2011.02.001,http://www.sciencedirect.com/ science/article/pii/S2210650211000022. [124] L.C.T. Bezerra, M. López-Ibáñez, T. Stützle, Automatic component-wise design of multiobjective evolutionary algorithms, IEEE Trans. Evol. Comput. 20 (3) (2016) 403–417, http://dx.doi.org/10.1109/TEVC.2015.2474158. [125] J. Derrac, S. García, D. Molina, F. Herrera, A practical tutorial on the use of nonparametric statistical tests as a methodology for comparing evolutionary and swarm intelligence algorithms, Swarm Evol. Comput. 1 (1) (2011) 3–18. [126] G. Vrbančič, L. Brezočnik, U. Mlakar, D. Fister, I. Fister Jr, NiaPy: Python Microframework for building nature-inspired algorithms, J. Open Source Softw. 3 (2018). [127] A. Benitez-Hidalgo, A.J. Nebro, J. Garcia-Nieto, I. Oregi, J. Del Ser, jMetalPy: A Python Framework for Multi-Objective Optimization with Metaheuristics, 2019, arXiv:1903.02915.