BACKGROUND
with the thriving development and confirmative significance of computational chemistry in drug discovery, more and more medicinal chemists and pharmacologists are using computational methods in their drug discovery research <cit> , and numerous drugs developed based on the clues provided by computations  have entered clinical trials or have been launched into the market already <cit> .
for the computational chemist, an attractive goal is to develop computer programs capable of automatically evaluating large-scale chemical libraries .
these computational methods are referred to as virtual screening  <cit> .
in general, two strategies have been employed in virtual screening: , using pharmacophore-based database searching  methods to identify potential hits from chemical libraries, mostly in the cases where three-dimensional  structures of the targets are unknown; and , using molecular docking approaches to screen the libraries in cases where the 3d structures of the targets are available <cit> .
because more and more 3d structures of drug target proteins are available, vs with molecular docking approaches has become promising, as demonstrated by numerous recent examples <cit> .
the core steps of structure-based virtual screening  are docking and scoring.
since kuntz et al.
<cit>  published the first docking algorithm dock in  <dig>  numerous docking programs have been developed during the past two decades  <cit> .
several comprehensive reviews of the advances of docking algorithms and applications have been published  <cit> .
scoring  the compounds retrieved from a database is performed simultaneously with the docking simulation.
molecular docking is a typical optimization problem, for it is difficult to obtain the global optimum solution.
as the fitness during the optimization process, scoring function should be fast and accurate enough, allowing simultaneous ranking of the retrieved poses in the optimization process.
based on this idea, several scoring functions have been developed  <cit> .
unfortunately, there is no scoring function developed so far that can reliably and consistently predict a ligand-protein binding mode and the binding affinity at the same time <cit> .
therefore, heuristic docking and consensus score strategies are often used in virtual screening  <cit> .
since a huge number of compounds in a database have to be thoroughly tested in the virtual screening process, several crucial issues have to be addressed.
one is the computational cost; only docking programs capable of docking a flexible ligand within a reasonable time scale are acceptable for virtual screening <cit> .
the other one is the ability to discriminate between true actives and inactive compounds; only those docking approaches able to distinguish the active molecules rapidly and accurately, are suitable for virtual screening applications in practice.
although the consensus score strategy has demonstrated an advantage over single scores, it is actually based on the results of single scoring optimization.
this means that consensus scoring only re-scores a limited molecules or conformations , thereby inevitably losing a number of true positives <cit> .
to some extent, consensus scoring seems to be far-fetched and artificial  <cit> .
most conformational optimization methods in docking program can only deal with a single objective, such as binding energy, shape complementarity, or chemical complementarity.
however, real-world problems normally involve multiple objectives  or optimization criteria, which should be satisfied simultaneously, and suitable solutions to the overall problem cannot be found by using individual optimization algorithms with single objectives <cit> .
for example, an optimization solution for the binding affinity  between a ligand and a receptor is usually not the optimization solution for other criteria .
similar problems in combinatorial library design and structural superposition of three-dimensional molecules have been reported <cit> .
thus, there is a need for an optimization algorithm compromising several objectives, which may result in more reasonable and robust binding modes between ligands and receptors.
in fact, it is a problem of multi-objective optimization  <cit> , which tends to find a set of alternative good compromises, generally known as pareto-optimal solutions.
these solutions are optimal because no other solutions in the search space superior to them when all objectives are considered.
then, the 'optimum' is chosen by the design which fits better in a certain application.
<cit>  there are more than twenty mathematical multi-objective optimization techniques <cit> .
however, due to their inherent parallelism, evolutionary algorithms  and genetic algorithms  are still the top priority in terms of finding multiple pareto-optimal solutions for multi-objective optimization problems <cit> .
in this paper, two sets of mo methods, denoted mosfom , were adopted for the binding conformation search of a small molecule within the binding site of a protein using two scoring functions as the objectives.
different from consensus score, mosfom does not re-score or re-rank the candidates from the primary virtual screening with one or several other scoring functions, but scores all the molecules in a chemical database with two or more scoring functions simultaneously during the binding pose optimization.
testing results indicate that mosfom, which is able to enhance hit rates and greatly reduces the false-positive rate, is more robust and reasonable than the consensus strategy as an alternate tool for large-scale library virtual screening.
here, mosfom emphasizes a new strategy to obtain the most reasonable binding conformation <cit>  and increase hit rates with several scoring functions rather than to accurately predict the binding free energy or the combination of several scoring functions.
consequently, mosfom can be used in the prioritization of ligands in high-throughput virtual screening.
methods
preparation of target proteins
thrombin, the estrogen receptor alpha, and cyclooxygenase- <dig>  have been used as target proteins for testing the newly developed docking algorithms in this study.
the coordinates of the x-ray crystal structures of these three proteins were retrieved from the protein data bank  <cit> , including thrombin in complex with mqpa at  <dig>  Å resolution  <cit> , estrogen receptor in complex with 4-hydroxytamoxifen at  <dig>  Å resolution  <cit> , and cyclooxygenase- <dig> in complex with sc- <dig> at  <dig>  Å resolution  <cit> .
all water molecules were removed from the protein structures.
after extraction of bound ligands, all hydrogen atoms and the kollman all-atom charges were assigned to the proteins using the biopolymer module of sybyl v <dig>  .
finally, for each protein target, the binding site was defined as the residues around the bound ligand within  <dig>  Å. gasteiger and marsili charges <cit>  were assigned to the extracted ligand of each protein.
preparation of compounds libraries
for all three targets , the active compounds were selected from the mddr .
active compounds with molecular weights between  <dig> and  <dig> were selected as drug-like compounds, and those containing water molecules and ions were eliminated.
table  <dig> shows the number of active compounds for each target.
another  <dig>  randomly 'varied' compounds with molecular weights between  <dig> and  <dig> were selected as drug-concerned decoys from the available chemical directory  after eliminating chemical reagents and inorganic compounds by means of corina.
all compounds of different test libraries were stored as sdf format file using the mdl isis_base program , and their three-dimensional coordinates were converted using a script written in the sybyl programming language , and gasteiger-marsili atomic charges were assigned to each molecule.
the final coordinates of each molecule were then stored in multi-mol <dig> files.
protonation states were correctly given for all the compounds.
scoring functions
the energy score and contact score of dock  <cit>  were used in this work.
energy score based on the amber force field <cit>  was composed of steric and electrostatic terms.
ε = 4r was used for the coefficient of the dielectric for the coulomb potential, and lennard-jones 6– <dig> was used for van der waals  potential.
contact score is a summation of the heavy atom contacts between the ligand and the receptor, which provides a simple assessment of shape complementarity; if two atoms approach close enough to bump into one another, then the interaction can be penalized by a certain amount.
in this study, a  <dig> -Å contact distance cutoff was used, with  <dig> penalized for each clash.
calculation of the enrichment factor
the enrichment factor is a key parameter to evaluate the quality of the docking and scoring compared to a random selection <cit> .
the enrichment factor  is defined as
ef=hitssns/hitstnt
where hitss is the number of active compounds in the sampled subset, hitst is the total number of active compounds in the database and n is the number of compounds.
in general, the enrichment factor is against random screening; it is evident that the maximum enrichment is determined by the total number of active compounds and the total number of molecules in the database.
for instance, there are  <dig> active compounds among the total  <dig> molecules in the database for the cox- <dig> case, i.e.
the achievable maximum is 10695/ <dig> ≈  <dig>  if 5% of active compounds were found among the top 1% of the database, then the enrichment factor would be fivefold over random  at the 1% of the database.
there are three criteria to evaluate the effectiveness of a docking program as indicated by wei et al.
<cit> : the value, the location of the vertex, and the percentage of active compounds retrieved, which represent the ability to find active compound of the docking program and scoring function.
ε-moea method
a steady-state moea based on the ε-dominance concept <cit>  is a pragmatic and fast multi-objective evolutionary algorithm for finding well-spread pareto-optimal solutions.
an ea population p and an archival population e  were used as two co-evolving populations in ε-moea.
after initialization, two solutions from p and e were chosen for mating.
then, each of these offspring solutions was compared with the archive and the ea population for possible inclusion.
for the case of j objectives, the search space and objective space were divided into hyper-boxes  to ensure that a hyper-box could be occupied by only one solution through comparing with an identification array of the archive, which guarantees the diversity of the archive.
the total number of pareto-optimal solutions, i.e.
the final size of suitable solutions, can be controlled approximately by adjusting the εj value .
two scoring functions, the energy score and the contact score of dock, were considered in practical virtual screening.
traditional optimization methods, such as the simplex method used in dock, are unsuitable for multi-objective optimization.
for the impartiality of the comparison, gasdock, which also uses a stochastic algorithm based on modified multi-population genetic algorithm <cit> , was adopted for single scoring function optimization.
the contact score also performed reasonably well in this study, which coincides with other studies  <cit> .
energy score and contact score were considered as the two objectives in ε-moea , that is, a set of pareto-optimal solutions, which were satisfied simultaneously with energy score and shape complementarity, were obtained consequently .
the multi-objective optimization model of ε-moea consists of a set of n parameters , a set of l objective functions, and a set of m constraints.
objective functions and constraints are functions of the decision variables.
it can be formulated mathematically as follows:
miny=f=,f <dig> ...,fl)s.t.g=,g <dig> ...,gm)≤0wherex=∈xy=∈y
and x is the design vector x = {tx, ty, tz, rx, ry, rz, tb <dig>  ⋯ tbn}t, in which  and  are respectively the state variables of translation and rotation of the entire ligand for the orientation search; and tb <dig>  ⋯, tbn are the torsion angles of the n rotatable bonds of the ligand for the conformation search.
accordingly, the constraints for the design variables s) can be represented as
{x¯≤tx≤x¯ y¯≤ty≤y¯ z¯≤tz≤z¯−π≤rx,ry,rz,tb <dig> ⋯,tbn≤π
y is the objective vector, which consists of energy score f <dig> and contact score f <dig>  respectively.
of course, more scoring function can be used in this method, but no more scoring function source can be accessed.
here the scoring functions of dock were just used to deal with this problem.
x is denoted as the decision space, and y is called the objective space.
selection of the optimum and ranking in ε-moea case
as stated above, a set of pareto-optimal solutions were obtained by using ε-moea, which simultaneously satisfied energy score and shape complementarity.
there are two ways to select an optimal solution from the set of pareto-optimal solutions: moea_nrg or 'energy score ≻ contact score' with the lowest energy conformation and acceptable shape complementarity; and moea_cnt or 'contact score ≻ energy score' with the best shape complementarity conformation and acceptable energy score.
the results of moea_nrg and moea_cnt were all compared with their corresponding individual scoring functions, respectively
efmoga
briefly, a new fast flexible docking program   <cit>  was developed using a multi-population genetic algorithm based on information entropy <cit> .
in comparison with other docking methods, information entropy was employed in the genetic algorithm of gasdock and the size of the narrowed space was used as the convergence criterion, ensuring that gasdock can converge rapidly and steadily.
a detailed description of the algorithm has been presented elsewhere <cit> .
in this paper, efmoga-based gasdock was applied to solve the above-mentioned multiple scoring function problem.
according to gasdock, the optimization problem  can be transformed into the following evaluation function model
minh=1sln⁡∑i=1lexp⁡)s.t.g=,g <dig> ...,gm)≤0wherex=∈xy=∈y
where λi is the weight of each objective.
).
by varying the weights, a set of noninferior solutions was generated.
s is a positive real variable, when s → ∞, the minimization problem converges to the maximal fi, and the minimum is the pareto solution
a knowledge-based method was adopted to obtain appropriate weights for virtual screening, i.e.
a set of conformations of the ligand, which are obtained through adjusting the weights, is compared with the ligand in the x-ray crystal structure, and the weights with minimum rmsd are selected.
these are considered to best reflect the test values for the current target system, i.e.
the most reasonable solution.
during this process, useful information can be acquired simultaneously such as the adaptability of different scoring functions to the current system, that is, which scoring function contributes more to the current system, and which scoring function is more sensitive.
as a result, various weights are obtained for different test systems in virtual screening .
for virtual screening of a database, each system only adopts one set of weights, which makes the comparison more equitable.
different from ε-moea, efmoga yields one solution corresponding to the weight in the pareto-optimal solutions rather than a set of pareto-optimal solutions.
certainly, compared with single objective optimization, the solution obtained here may be worse than any single objective solution, but more reasonable than any extremum obtained by a single objective method because it has considered multiple-objective functions.
during the preparation of this manuscript, grosdidier et al.
have developed a new docking software, eadock, based on a multi-objective optimization algorithm, and the results indicate that eadock can accurately predict binding modes for ligand-protein complexes with two fitness functions <cit> .
RESULTS
as mentioned above, our study is not aimed at calculating absolute values for the free energy of binding and for the affinity, but focuses on the ranking of acceptable conformations for the multiple solution space.
because of the different preferences in the selection of pareto-optimal solutions with ε-moea, the results of three different approaches will be compared with individual scoring functions, namely one with preferred energy score for ε-moea , the other with preferred contact score for ε-moea , and the third for efmoga.
docking to the buried lipophilic site of cox-2
the cox- <dig> ligand binding site is a completely buried, narrow, confined and predominately lipophilic cavity  <cit> , which does not accommodate many orientations or conformations of a ligand.
consequently, shape complementarity was expected to be very important.
moea_nrg efficiently eliminates conformations that could not satisfy shape complementarity although they displayed good energy scores.
the maximum enrichment was  <dig> -fold over random, which was reached at the top  <dig> % of the database when using a single energy score for optimization, but the maximum enrichment was  <dig> -fold over random, and moved up to the top  <dig> % of the database when using moea_nrg .
moea_cnt also got a good enrichment  at the top  <dig> % of the database, but it was worse than single contact score with ef =  <dig>  at the top  <dig> %.
there are several possible explanations; first, the conformation with highest contact score cannot satisfy energy score, that is, the conformation that has the best contact score but a bad energy score not among the pareto-optimal solutions will be eliminated; second, the scoring function is imprecise which is the common disadvantage for all the scoring functions.
interestingly, although contact score maintained higher enrichment among the top 2% of the database, the same enrichment as that of moea_nrg was obtained at the top 2% of the database, and from this point, the enrichment peak of moea_nrg became broader.
efmoga has a similar curve as that of moea_cnt, and a maximum of  <dig> -fold over random was reached at the top  <dig> % of the database.
in the real application, only a small number  was our interesting section within a large source library.
moea_nrg exhibited a good performance in this case , moea_nrg maintained the best performance among all the methods, and retrieved  <dig> % of the cox- <dig> active compounds within the top 2% of the total library .
efmoga outperformed single energy score among the top 5% of the library, and retrieved 6% of active compounds among the top 2% of the database, although it performed worse than single contact score.
whereas, moea_cnt represents a relatively weak ability than single contact score before the top 5% of the database, but greatly increases thereafter.
docking to the relatively large hydrophobic site of estrogen receptor
as an example of a well-understood target with therapeutic relevance whose crystallographic structural data were publicly available, estrogen receptor was chosen.
like cox- <dig>  the binding site of estrogen receptor is a relatively larger, fully-enclosed lipophilic cavity , which is little opened to solvent, there are acceptor groups at either end that can form hydrogen bonds with ligands, but it is predominantly hydrophobic on the whole <cit> .
for this case, efmoga performed well among the top 5% of the database, especially before the top 2%, and reached the maximum enrichment  at the top  <dig> % of the database .
single energy score exhibited the poorest performance in all methods, the maximum enrichment is less than fivefold over random at the top  <dig> % of the database, although moea_nrg reached the same enrichment, moea_nrg outperformed single energy score after that.
different from cox- <dig>  moea_cnt indicated a very good performance than single contact score for this case, it is more reasonable that we preferred contact score from pareto-optimal solutions to select the better shape satisfaction with slight difference in energy score.
using single contact score, the maximum enrichment  was reached at the top  <dig> % of the database, but the maximum enrichment rose to  <dig> , and moved up to the top  <dig> % of the database when using moea_cnt.
similar to cox- <dig>  approximately  <dig> % of active compounds were retrieved using moea_nrg at the top 2% of the database , which was about twofold over that of single energy score.
relative to single contact score, moea_cnt retrieved  <dig> % of the active compounds at the top 2% of the database, but only less than 5% of the active compounds were retrieved using single contact score.
it is encouraging that efmoga represents better performance than energy score and contact score at the top 2% of the database, with 8% of active compounds retrieved.
docking to the intermediate polarity site of thrombin
in contrast to cox- <dig> and estrogen receptor, the binding site of thrombin is an intermediate polarity site <cit> , with a large hydrophobic pocket  which is a more exposed binding site to solvent <cit> , in addition, there is a s <dig> specificity pocket, which is a narrow and restricted pocket comprising the carboxylate group of asp <dig> at the bottom, most thrombin inhibitors form hydrogen bonds with asp <dig> and also to gly <dig> .
due to many polar groups in the binding site, it is difficult to distinguish active compounds for dock energy score which is most reliable for the apolar active site, unlike those scoring functions with extra consideration of hydrogen bonding interactions such as pmf or flexx score <cit> .
as expected, although the maximum enrichment was  <dig> at the top  <dig> % of the database for moea_nrg, and the maximum enrichment was  <dig>  at the top  <dig> % of the database when using single energy score, there is a slight increase in quantity using moea_nrg against single energy score at the top 2% of the database .
moea_nrg slightly underperformed than single energy score after the top 2% of the database .
like er case, moea_cnt performed the best in all the methods, the maximum enrichment  was reached at the top  <dig> % of the database, approximately twofold over contact score, with a maximum value of  <dig>  at the top  <dig> % of the database , at the same time,  <dig> % of the active compounds were retrieved with moea_cnt at the top 2% of the database .
same as other two systems, efmoga continues a good performance at the top 2% of the database, with a  <dig>  enrichment reached at the top  <dig> % of the database, and retrieved approximately 5% of the active compounds at the top 2% of the database, which is less than that of contact score, but better than that of energy score, efmoga and moea_nrg performed not so well compared with single score function after the top 5% of the database.
discussion
performance of mosfom and characteristics of binding site and preferential selection of scoring function
simple as it is, contact score performed unexpectedly well for all the three systems in our study.
especially at the top 5%, contact score performed better than energy score, possibly arising from molecular diversity of the randomly selected acd database, which distributes well not only in heavy atoms but also molecular torsions.
contact score can rapidly seek out those whose geometry shape can satisfy the binding site , which is especially obvious to those with completely buried and narrow cavity.
however, energy score presented a weak ability to distinguish the active compounds from the decoys with shape satisfaction.
it is encouraging that compared with individual score, every method of mosfom presented an increase not only in enrichment but also hit rates of the active compounds retrieved in some ways.
efmoga with different weights in different systems obtained better results than single score at the top 5%  of the database for all three cases in this study , but gradually faded afterwards.
it can be easily understood that those active compounds with high affinity and good shape complementarity were first retrieved, and then those emphasizing particularly on the larger weight in different test system were obtained, consecutively.
therefore, if not knowing which scoring function will work better for the test system in advance, more enrichment can be obtained with efmoga in the most interesting section  of the database.
it should be noted that efmoga, with different weights for different test systems through experimental protein-ligand complex, is more reasonable and different from the simple linear combination  <cit>  of several scoring functions or screenscore through a combination of scoring terms <cit> .
efmoga seems more pertinent and advantageous, and will be a tendency for scoring and ranking to consider different scoring functions focusing on diverse aspects for different systems.
different from efmoga, ε-moea selected the preferential solution from a pareto-optimal solutions family  through multi-objective optimization using multiple scoring functions .
with different preferences, moea_nrg and moea_cnt all outperformed their corresponding individual scoring function, i.e., moea_nrg corresponding to energy score and moea_cnt corresponding to contact score.
a higher enrichment and hit rates of the active compounds were obtained using moea_nrg than using individual energy score in cox- <dig> and er system, but there is an inverse phenomenon for thrombin case.
there are several strong hydrogen bond interactions between the ligand and the residues in s <dig> pocket for thrombin system.
if we select the one with lowest energy score from the docked conformations with acceptable shape complementarity, more compounds will have good shape complementarity conformations but are only slightly different in energy score, since the energy score does not pay more additional attention to the h-bond interaction.
in this case, it will be more difficult to distinguish the active compounds.
on the contrary, moea_cnt performed well for the thrombin and er case, where an open cavity exists, accommodating enough orientations or conformations of a ligand, therefore, it seems more reasonable that moea_cnt selected the best shape complementarity conformation with differences in energy score.
so, if one does know which scoring function will work better for the test system in advance, more enrichment and better performance can be obtained using ε-moea through relative preference.
limitations of mosfom, further improvements and prospect in bioinformatics and drug design
more scoring functions, which pay special attention to hydrogen bonding interaction, or chemical complementarity, solvent effect, should be involved in multi-objective optimization, however, there are only energy score and contact score available for us.
although mosfom performed better than individual score, there will be more choices for mofsom in specific test system if more interactions were considered.
we can determine different weights of scoring functions because efmoga can be determined to deal with different cases and maybe prefer others  for ε-moea if we know which scoring function will work well for specific test system in advance.
therefore, to develop scoring functions with focus on hydrogen bond, icon or solvent effect and to utilize multi-objective optimization to acquire more reasonable conformations is one of our future research interest.
at the same time, ignorance of the impact of protein flexibility is another limitation for mosfom, small changes of the receptor flexibility result in larger variety of binding affinities <cit> , and docking to a single receptor conformation will significantly reduce the chance of finding the correct pose.
considering protein flexibility to acquire correct protein-ligand binding conformation is crucial for medicinal chemists to find out the structure-activity relationship.
we will fulfil this work in the future.
also, computational time of multi-objective optimization is another issue to be solved.
there are numerous multi-objective optimization methods, however, they are not suitable for large-scale virtual screening because of the huge time consumption, in this paper, ε-moea, which is a steady-state moea based on the ε-dominance concept, fulfilled virtual screening quickly within about one minute for one molecule, and efmoga, which is a multi-population moga based on information entropy, accomplished one docking within  <dig> seconds using one cpu on sgi origin <dig> hardware.
however, more improvements are needed to find a set of well-distributed solutions close to the true pareto frontier pragmatically and efficiently.
as stated above, multi-objective optimization approach has been adopted in the library design program moselect <cit>  and quantitative structure activity relationships program moqsar <cit> .
we believe that moea or moga also can be applied to biological sequence alignment, protein fold recognition, conformational generation and adme/tox  with more and more factors taken into account in bioinformatics and drug design.
CONCLUSIONS
development of a fast and accurate scoring function in virtual screening is still a hot and pending issue in current research, different scoring functions focus on different aspects of ligand binding, and no single scoring can satisfy all the systems well, therefore, consensus score was put forward.
<cit>  with several other scoring functions, consensus score re-assessed those conformations optimized using a primary scoring function, but it is not really robust from the viewpoint of optimization.
all of these give rise to another heuristic thinking for us, is it possible and more rational to find a most reasonable conformation in the process of optimization using two or more scoring functions in virtual screening?
since present scoring functions can not satisfy all the cases, multi-objective optimization method can be a good compromise.
in this study, we present two kinds of multi-objective optimization method, called mosfom, in virtual screening which simultaneously considers energy score and contact score.
a set of pareto-optimal solutions were obtained which can simultaneously satisfy energy and shape complementarity using ε-moea method, then through two different preferences, the binding conformations of pareto-optimal solutions with lowest energy scoring or best shape complementarity were ranked among all selected conformations.
however, efmoga acquires different weights of each scoring function for different systems based on experimental x-ray complex structure, that is, for different system, efmoga regards the varying weight as the degree of contribution of individual scoring function, which generally focuses on one aspect for a special type of binding sites, this means the scoring function with bigger weight is more suitable for the binding site.
we used  <dig> random compounds as the decoys, active compounds selected from mddr database were randomly merged into the decoys, respectively.
to ensure justness of comparison, gasdock, based on an improved genetic algorithm, was used as a benchmark for single-objective optimization.
results show that mosfom can enhance the enrichment and greatly increase the hit rates compared with individual score .
for three different kinds of binding sites, mosfom represents excellent ability of distinction of active compounds with energy and shape complementarity.
efmoga specially performed well at the top 2% of the database, moea_nrg and moea_cnt performed better than individual scoring function if a proper type of binding site was selected.
multi-objective optimization method was successfully applied in virtual screening with two different scoring functions, which can gain reasonable binding pose, and can be furthermore ranked with those potentially compromised conformations of each compound which abandon those conformations that can not satisfy overall objective functions.
by analyzing the characteristics of binding site in advance, the most effective multi-objective optimization method is adopted, which is meaningful for all current scoring functions since they can not suit all cases.
more specific scoring functions and protein flexibility will be taken into account in our future work.
authors' contributions
hl, hz and jl contributed to the development and validation of the method, mz, lk, and xl contributed to the analysis and data interpretation, xw and hj conceived the idea of the method, and provided guidance for its development and revised the subsequent drafts of this manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
docking example and method proof.
a molecular docking example for reproduction of thymidine kinase  complex, and the mathematical proof for evaluation function multi-objective optimization genetic algorithm .
click here for file
acknowledgements
the authors thank prof. k. deb for the source code of ε-moea.
we also thank prof. rolf hilgenfeld for careful reading this manuscript, and the referees for helpful suggestions and comments.
this work was supported by the special fund for major state basic research project , the national natural science foundation of china , the shanghai committee of science and technology, and the  <dig> hi-tech program of china .
h.l.
was also supported by knowledge innovation program of the chinese academy of sciences .

 with the thriving development and confirmative significance of computational chemistry in drug discovery, more and more medicinal chemists and pharmacologists are using computational methods in their drug discovery research, and numerous drugs developed based on the clues provided by computations have entered clinical trials or have been launched into the market already. for the computational chemist, an attractive goal is to develop computer programs capable of automatically evaluating large-scale chemical libraries. these computational methods are referred to as virtual screening. in general, two strategies have been employed in virtual screening:, using pharmacophore-based database searching methods to identify potential hits from chemical libraries, mostly in the cases where three-dimensional structures of the targets are unknown; and, using molecular docking approaches to screen the libraries in cases where the 3d structures of the targets are available. because more and more 3d structures of drug target proteins are available, vs with molecular docking approaches has become promising, as demonstrated by numerous recent examples. the core steps of structure-based virtual screening are docking and scoring. since kuntz et al. published the first docking algorithm dock in numerous docking programs have been developed during the past two decades. several comprehensive reviews of the advances of docking algorithms and applications have been published. scoring the compounds retrieved from a database is performed simultaneously with the docking simulation. molecular docking is a typical optimization problem, for it is difficult to obtain the global optimum solution. as the fitness during the optimization process, scoring function should be fast and accurate enough, allowing simultaneous ranking of the retrieved poses in the optimization process. based on this idea, several scoring functions have been developed. unfortunately, there is no scoring function developed so far that can reliably and consistently predict a ligand-protein binding mode and the binding affinity at the same time. therefore, heuristic docking and consensus score strategies are often used in virtual screening. since a huge number of compounds in a database have to be thoroughly tested in the virtual screening process, several crucial issues have to be addressed. one is the computational cost; only docking programs capable of docking a flexible ligand within a reasonable time scale are acceptable for virtual screening. the other one is the ability to discriminate between true actives and inactive compounds; only those docking approaches able to distinguish the active molecules rapidly and accurately, are suitable for virtual screening applications in practice. although the consensus score strategy has demonstrated an advantage over single scores, it is actually based on the results of single scoring optimization. this means that consensus scoring only re-scores a limited molecules or conformations, thereby inevitably losing a number of true positives. to some extent, consensus scoring seems to be far-fetched and artificial. most conformational optimization methods in docking program can only deal with a single objective, such as binding energy, shape complementarity, or chemical complementarity. however, real-world problems normally involve multiple objectives or optimization criteria, which should be satisfied simultaneously, and suitable solutions to the overall problem cannot be found by using individual optimization algorithms with single objectives. for example, an optimization solution for the binding affinity between a ligand and a receptor is usually not the optimization solution for other criteria. similar problems in combinatorial library design and structural superposition of three-dimensional molecules have been reported. thus, there is a need for an optimization algorithm compromising several objectives, which may result in more reasonable and robust binding modes between ligands and receptors. in fact, it is a problem of multi-objective optimization, which tends to find a set of alternative good compromises, generally known as pareto-optimal solutions. these solutions are optimal because no other solutions in the search space superior to them when all objectives are considered. then, the'optimum' is chosen by the design which fits better in a certain application. there are more than twenty mathematical multi-objective optimization techniques. however, due to their inherent parallelism, evolutionary algorithms and genetic algorithms are still the top priority in terms of finding multiple pareto-optimal solutions for multi-objective optimization problems. in this paper, two sets of mo methods, denoted mosfom, were adopted for the binding conformation search of a small molecule within the binding site of a protein using two scoring functions as the objectives. different from consensus score, mosfom does not re-score or re-rank the candidates from the primary virtual screening with one or several other scoring functions, but scores all the molecules in a chemical database with two or more scoring functions simultaneously during the binding pose optimization. testing results indicate that mosfom, which is able to enhance hit rates and greatly reduces the false-positive rate, is more robust and reasonable than the consensus strategy as an alternate tool for large-scale library virtual screening. here, mosfom emphasizes a new strategy to obtain the most reasonable binding conformation and increase hit rates with several scoring functions rather than to accurately predict the binding free energy or the combination of several scoring functions. consequently, mosfom can be used in the prioritization of ligands in high-throughput virtual screening. methods preparation of target proteins thrombin, the estrogen receptor alpha, and cyclooxygenase- have been used as target proteins for testing the newly developed docking algorithms in this study. the coordinates of the x-ray crystal structures of these three proteins were retrieved from the protein data bank, including thrombin in complex with mqpa at Å resolution, estrogen receptor in complex with 4-hydroxytamoxifen at Å resolution, and cyclooxygenase- in complex with sc- at Å resolution. all water molecules were removed from the protein structures. after extraction of bound ligands, all hydrogen atoms and the kollman all-atom charges were assigned to the proteins using the biopolymer module of sybyl v. finally, for each protein target, the binding site was defined as the residues around the bound ligand within Å. gasteiger and marsili charges were assigned to the extracted ligand of each protein. preparation of compounds libraries for all three targets, the active compounds were selected from the mddr. active compounds with molecular weights between and were selected as drug-like compounds, and those containing water molecules and ions were eliminated. table shows the number of active compounds for each target. another randomly'varied' compounds with molecular weights between and were selected as drug-concerned decoys from the available chemical directory after eliminating chemical reagents and inorganic compounds by means of corina. all compounds of different test libraries were stored as sdf format file using the mdl isis_base program, and their three-dimensional coordinates were converted using a script written in the sybyl programming language, and gasteiger-marsili atomic charges were assigned to each molecule. the final coordinates of each molecule were then stored in multi-mol files. protonation states were correctly given for all the compounds. scoring functions the energy score and contact score of dock were used in this work. energy score based on the amber force field was composed of steric and electrostatic terms. ε = 4r was used for the coefficient of the dielectric for the coulomb potential, and lennard-jones 6– was used for van der waals potential. contact score is a summation of the heavy atom contacts between the ligand and the receptor, which provides a simple assessment of shape complementarity; if two atoms approach close enough to bump into one another, then the interaction can be penalized by a certain amount. in this study, a -Å contact distance cutoff was used, with penalized for each clash. calculation of the enrichment factor the enrichment factor is a key parameter to evaluate the quality of the docking and scoring compared to a random selection. the enrichment factor is defined as ef=hitssns/hitstnt where hitss is the number of active compounds in the sampled subset, hitst is the total number of active compounds in the database and n is the number of compounds. in general, the enrichment factor is against random screening; it is evident that the maximum enrichment is determined by the total number of active compounds and the total number of molecules in the database. for instance, there are active compounds among the total molecules in the database for the cox- case, i.e. the achievable maximum is 10695/ ≈ if 5% of active compounds were found among the top 1% of the database, then the enrichment factor would be fivefold over random at the 1% of the database. there are three criteria to evaluate the effectiveness of a docking program as indicated by wei et al.: the value, the location of the vertex, and the percentage of active compounds retrieved, which represent the ability to find active compound of the docking program and scoring function. ε-moea method a steady-state moea based on the ε-dominance concept is a pragmatic and fast multi-objective evolutionary algorithm for finding well-spread pareto-optimal solutions. an ea population p and an archival population e were used as two co-evolving populations in ε-moea. after initialization, two solutions from p and e were chosen for mating. then, each of these offspring solutions was compared with the archive and the ea population for possible inclusion. for the case of j objectives, the search space and objective space were divided into hyper-boxes to ensure that a hyper-box could be occupied by only one solution through comparing with an identification array of the archive, which guarantees the diversity of the archive. the total number of pareto-optimal solutions, i.e. the final size of suitable solutions, can be controlled approximately by adjusting the εj value. two scoring functions, the energy score and the contact score of dock, were considered in practical virtual screening. traditional optimization methods, such as the simplex method used in dock, are unsuitable for multi-objective optimization. for the impartiality of the comparison, gasdock, which also uses a stochastic algorithm based on modified multi-population genetic algorithm, was adopted for single scoring function optimization. the contact score also performed reasonably well in this study, which coincides with other studies. energy score and contact score were considered as the two objectives in ε-moea, that is, a set of pareto-optimal solutions, which were satisfied simultaneously with energy score and shape complementarity, were obtained consequently. the multi-objective optimization model of ε-moea consists of a set of n parameters, a set of l objective functions, and a set of m constraints. objective functions and constraints are functions of the decision variables. it can be formulated mathematically as follows: miny=f=,f...,fl)s.t.g=,g...,gm)≤0wherex=∈xy=∈y and x is the design vector x = {tx, ty, tz, rx, ry, rz, tb ⋯ tbn}t, in which and are respectively the state variables of translation and rotation of the entire ligand for the orientation search; and tb ⋯, tbn are the torsion angles of the n rotatable bonds of the ligand for the conformation search. accordingly, the constraints for the design variables s) can be represented as {x¯≤tx≤x¯ y¯≤ty≤y¯ z¯≤tz≤z¯−π≤rx,ry,rz,tb ⋯,tbn≤π y is the objective vector, which consists of energy score f and contact score f respectively. of course, more scoring function can be used in this method, but no more scoring function source can be accessed. here the scoring functions of dock were just used to deal with this problem. x is denoted as the decision space, and y is called the objective space. selection of the optimum and ranking in ε-moea case as stated above, a set of pareto-optimal solutions were obtained by using ε-moea, which simultaneously satisfied energy score and shape complementarity. there are two ways to select an optimal solution from the set of pareto-optimal solutions: moea_nrg or'energy score ≻ contact score' with the lowest energy conformation and acceptable shape complementarity; and moea_cnt or'contact score ≻ energy score' with the best shape complementarity conformation and acceptable energy score. the results of moea_nrg and moea_cnt were all compared with their corresponding individual scoring functions, respectively efmoga briefly, a new fast flexible docking program was developed using a multi-population genetic algorithm based on information entropy. in comparison with other docking methods, information entropy was employed in the genetic algorithm of gasdock and the size of the narrowed space was used as the convergence criterion, ensuring that gasdock can converge rapidly and steadily. a detailed description of the algorithm has been presented elsewhere. in this paper, efmoga-based gasdock was applied to solve the above-mentioned multiple scoring function problem. according to gasdock, the optimization problem can be transformed into the following evaluation function model minh=1sln⁡∑i=1lexp⁡)s.t.g=,g...,gm)≤0wherex=∈xy=∈y where λi is the weight of each objective. ). by varying the weights, a set of noninferior solutions was generated. s is a positive real variable, when s → ∞, the minimization problem converges to the maximal fi, and the minimum is the pareto solution a knowledge-based method was adopted to obtain appropriate weights for virtual screening, i.e. a set of conformations of the ligand, which are obtained through adjusting the weights, is compared with the ligand in the x-ray crystal structure, and the weights with minimum rmsd are selected. these are considered to best reflect the test values for the current target system, i.e. the most reasonable solution. during this process, useful information can be acquired simultaneously such as the adaptability of different scoring functions to the current system, that is, which scoring function contributes more to the current system, and which scoring function is more sensitive. as a result, various weights are obtained for different test systems in virtual screening. for virtual screening of a database, each system only adopts one set of weights, which makes the comparison more equitable. different from ε-moea, efmoga yields one solution corresponding to the weight in the pareto-optimal solutions rather than a set of pareto-optimal solutions. certainly, compared with single objective optimization, the solution obtained here may be worse than any single objective solution, but more reasonable than any extremum obtained by a single objective method because it has considered multiple-objective functions. during the preparation of this manuscript, grosdidier et al. have developed a new docking software, eadock, based on a multi-objective optimization algorithm, and the results indicate that eadock can accurately predict binding modes for ligand-protein complexes with two fitness functions.  as mentioned above, our study is not aimed at calculating absolute values for the free energy of binding and for the affinity, but focuses on the ranking of acceptable conformations for the multiple solution space. because of the different preferences in the selection of pareto-optimal solutions with ε-moea, the results of three different approaches will be compared with individual scoring functions, namely one with preferred energy score for ε-moea, the other with preferred contact score for ε-moea, and the third for efmoga. docking to the buried lipophilic site of cox-2 the cox- ligand binding site is a completely buried, narrow, confined and predominately lipophilic cavity, which does not accommodate many orientations or conformations of a ligand. consequently, shape complementarity was expected to be very important. moea_nrg efficiently eliminates conformations that could not satisfy shape complementarity although they displayed good energy scores. the maximum enrichment was -fold over random, which was reached at the top % of the database when using a single energy score for optimization, but the maximum enrichment was -fold over random, and moved up to the top % of the database when using moea_nrg. moea_cnt also got a good enrichment at the top % of the database, but it was worse than single contact score with ef = at the top %. there are several possible explanations; first, the conformation with highest contact score cannot satisfy energy score, that is, the conformation that has the best contact score but a bad energy score not among the pareto-optimal solutions will be eliminated; second, the scoring function is imprecise which is the common disadvantage for all the scoring functions. interestingly, although contact score maintained higher enrichment among the top 2% of the database, the same enrichment as that of moea_nrg was obtained at the top 2% of the database, and from this point, the enrichment peak of moea_nrg became broader. efmoga has a similar curve as that of moea_cnt, and a maximum of -fold over random was reached at the top % of the database. in the real application, only a small number was our interesting section within a large source library. moea_nrg exhibited a good performance in this case, moea_nrg maintained the best performance among all the methods, and retrieved % of the cox- active compounds within the top 2% of the total library. efmoga outperformed single energy score among the top 5% of the library, and retrieved 6% of active compounds among the top 2% of the database, although it performed worse than single contact score. whereas, moea_cnt represents a relatively weak ability than single contact score before the top 5% of the database, but greatly increases thereafter. docking to the relatively large hydrophobic site of estrogen receptor as an example of a well-understood target with therapeutic relevance whose crystallographic structural data were publicly available, estrogen receptor was chosen. like cox- the binding site of estrogen receptor is a relatively larger, fully-enclosed lipophilic cavity, which is little opened to solvent, there are acceptor groups at either end that can form hydrogen bonds with ligands, but it is predominantly hydrophobic on the whole. for this case, efmoga performed well among the top 5% of the database, especially before the top 2%, and reached the maximum enrichment at the top % of the database. single energy score exhibited the poorest performance in all methods, the maximum enrichment is less than fivefold over random at the top % of the database, although moea_nrg reached the same enrichment, moea_nrg outperformed single energy score after that. different from cox- moea_cnt indicated a very good performance than single contact score for this case, it is more reasonable that we preferred contact score from pareto-optimal solutions to select the better shape satisfaction with slight difference in energy score. using single contact score, the maximum enrichment was reached at the top % of the database, but the maximum enrichment rose to, and moved up to the top % of the database when using moea_cnt. similar to cox- approximately % of active compounds were retrieved using moea_nrg at the top 2% of the database, which was about twofold over that of single energy score. relative to single contact score, moea_cnt retrieved % of the active compounds at the top 2% of the database, but only less than 5% of the active compounds were retrieved using single contact score. it is encouraging that efmoga represents better performance than energy score and contact score at the top 2% of the database, with 8% of active compounds retrieved. docking to the intermediate polarity site of thrombin in contrast to cox- and estrogen receptor, the binding site of thrombin is an intermediate polarity site, with a large hydrophobic pocket which is a more exposed binding site to solvent, in addition, there is a s specificity pocket, which is a narrow and restricted pocket comprising the carboxylate group of asp at the bottom, most thrombin inhibitors form hydrogen bonds with asp and also to gly. due to many polar groups in the binding site, it is difficult to distinguish active compounds for dock energy score which is most reliable for the apolar active site, unlike those scoring functions with extra consideration of hydrogen bonding interactions such as pmf or flexx score. as expected, although the maximum enrichment was at the top % of the database for moea_nrg, and the maximum enrichment was at the top % of the database when using single energy score, there is a slight increase in quantity using moea_nrg against single energy score at the top 2% of the database. moea_nrg slightly underperformed than single energy score after the top 2% of the database. like er case, moea_cnt performed the best in all the methods, the maximum enrichment was reached at the top % of the database, approximately twofold over contact score, with a maximum value of at the top % of the database, at the same time, % of the active compounds were retrieved with moea_cnt at the top 2% of the database. same as other two systems, efmoga continues a good performance at the top 2% of the database, with a enrichment reached at the top % of the database, and retrieved approximately 5% of the active compounds at the top 2% of the database, which is less than that of contact score, but better than that of energy score, efmoga and moea_nrg performed not so well compared with single score function after the top 5% of the database. discussion performance of mosfom and characteristics of binding site and preferential selection of scoring function simple as it is, contact score performed unexpectedly well for all the three systems in our study. especially at the top 5%, contact score performed better than energy score, possibly arising from molecular diversity of the randomly selected acd database, which distributes well not only in heavy atoms but also molecular torsions. contact score can rapidly seek out those whose geometry shape can satisfy the binding site, which is especially obvious to those with completely buried and narrow cavity. however, energy score presented a weak ability to distinguish the active compounds from the decoys with shape satisfaction. it is encouraging that compared with individual score, every method of mosfom presented an increase not only in enrichment but also hit rates of the active compounds retrieved in some ways. efmoga with different weights in different systems obtained better results than single score at the top 5% of the database for all three cases in this study, but gradually faded afterwards. it can be easily understood that those active compounds with high affinity and good shape complementarity were first retrieved, and then those emphasizing particularly on the larger weight in different test system were obtained, consecutively. therefore, if not knowing which scoring function will work better for the test system in advance, more enrichment can be obtained with efmoga in the most interesting section of the database. it should be noted that efmoga, with different weights for different test systems through experimental protein-ligand complex, is more reasonable and different from the simple linear combination of several scoring functions or screenscore through a combination of scoring terms. efmoga seems more pertinent and advantageous, and will be a tendency for scoring and ranking to consider different scoring functions focusing on diverse aspects for different systems. different from efmoga, ε-moea selected the preferential solution from a pareto-optimal solutions family through multi-objective optimization using multiple scoring functions. with different preferences, moea_nrg and moea_cnt all outperformed their corresponding individual scoring function, i.e., moea_nrg corresponding to energy score and moea_cnt corresponding to contact score. a higher enrichment and hit rates of the active compounds were obtained using moea_nrg than using individual energy score in cox- and er system, but there is an inverse phenomenon for thrombin case. there are several strong hydrogen bond interactions between the ligand and the residues in s pocket for thrombin system. if we select the one with lowest energy score from the docked conformations with acceptable shape complementarity, more compounds will have good shape complementarity conformations but are only slightly different in energy score, since the energy score does not pay more additional attention to the h-bond interaction. in this case, it will be more difficult to distinguish the active compounds. on the contrary, moea_cnt performed well for the thrombin and er case, where an open cavity exists, accommodating enough orientations or conformations of a ligand, therefore, it seems more reasonable that moea_cnt selected the best shape complementarity conformation with differences in energy score. so, if one does know which scoring function will work better for the test system in advance, more enrichment and better performance can be obtained using ε-moea through relative preference. limitations of mosfom, further improvements and prospect in bioinformatics and drug design more scoring functions, which pay special attention to hydrogen bonding interaction, or chemical complementarity, solvent effect, should be involved in multi-objective optimization, however, there are only energy score and contact score available for us. although mosfom performed better than individual score, there will be more choices for mofsom in specific test system if more interactions were considered. we can determine different weights of scoring functions because efmoga can be determined to deal with different cases and maybe prefer others for ε-moea if we know which scoring function will work well for specific test system in advance. therefore, to develop scoring functions with focus on hydrogen bond, icon or solvent effect and to utilize multi-objective optimization to acquire more reasonable conformations is one of our future research interest. at the same time, ignorance of the impact of protein flexibility is another limitation for mosfom, small changes of the receptor flexibility result in larger variety of binding affinities, and docking to a single receptor conformation will significantly reduce the chance of finding the correct pose. considering protein flexibility to acquire correct protein-ligand binding conformation is crucial for medicinal chemists to find out the structure-activity relationship. we will fulfil this work in the future. also, computational time of multi-objective optimization is another issue to be solved. there are numerous multi-objective optimization methods, however, they are not suitable for large-scale virtual screening because of the huge time consumption, in this paper, ε-moea, which is a steady-state moea based on the ε-dominance concept, fulfilled virtual screening quickly within about one minute for one molecule, and efmoga, which is a multi-population moga based on information entropy, accomplished one docking within seconds using one cpu on sgi origin hardware. however, more improvements are needed to find a set of well-distributed solutions close to the true pareto frontier pragmatically and efficiently. as stated above, multi-objective optimization approach has been adopted in the library design program moselect and quantitative structure activity relationships program moqsar. we believe that moea or moga also can be applied to biological sequence alignment, protein fold recognition, conformational generation and adme/tox with more and more factors taken into account in bioinformatics and drug design.  development of a fast and accurate scoring function in virtual screening is still a hot and pending issue in current research, different scoring functions focus on different aspects of ligand binding, and no single scoring can satisfy all the systems well, therefore, consensus score was put forward. with several other scoring functions, consensus score re-assessed those conformations optimized using a primary scoring function, but it is not really robust from the viewpoint of optimization. all of these give rise to another heuristic thinking for us, is it possible and more rational to find a most reasonable conformation in the process of optimization using two or more scoring functions in virtual screening? since present scoring functions can not satisfy all the cases, multi-objective optimization method can be a good compromise. in this study, we present two kinds of multi-objective optimization method, called mosfom, in virtual screening which simultaneously considers energy score and contact score. a set of pareto-optimal solutions were obtained which can simultaneously satisfy energy and shape complementarity using ε-moea method, then through two different preferences, the binding conformations of pareto-optimal solutions with lowest energy scoring or best shape complementarity were ranked among all selected conformations. however, efmoga acquires different weights of each scoring function for different systems based on experimental x-ray complex structure, that is, for different system, efmoga regards the varying weight as the degree of contribution of individual scoring function, which generally focuses on one aspect for a special type of binding sites, this means the scoring function with bigger weight is more suitable for the binding site. we used random compounds as the decoys, active compounds selected from mddr database were randomly merged into the decoys, respectively. to ensure justness of comparison, gasdock, based on an improved genetic algorithm, was used as a benchmark for single-objective optimization. results show that mosfom can enhance the enrichment and greatly increase the hit rates compared with individual score. for three different kinds of binding sites, mosfom represents excellent ability of distinction of active compounds with energy and shape complementarity. efmoga specially performed well at the top 2% of the database, moea_nrg and moea_cnt performed better than individual scoring function if a proper type of binding site was selected. multi-objective optimization method was successfully applied in virtual screening with two different scoring functions, which can gain reasonable binding pose, and can be furthermore ranked with those potentially compromised conformations of each compound which abandon those conformations that can not satisfy overall objective functions. by analyzing the characteristics of binding site in advance, the most effective multi-objective optimization method is adopted, which is meaningful for all current scoring functions since they can not suit all cases. more specific scoring functions and protein flexibility will be taken into account in our future work. authors' contributions hl, hz and jl contributed to the development and validation of the method, mz, lk, and xl contributed to the analysis and data interpretation, xw and hj conceived the idea of the method, and provided guidance for its development and revised the subsequent drafts of this manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 docking example and method proof. a molecular docking example for reproduction of thymidine kinase complex, and the mathematical proof for evaluation function multi-objective optimization genetic algorithm. click here for file acknowledgements the authors thank prof. k. deb for the source code of ε-moea. we also thank prof. rolf hilgenfeld for careful reading this manuscript, and the referees for helpful suggestions and comments. this work was supported by the special fund for major state basic research project, the national natural science foundation of china, the shanghai committee of science and technology, and the hi-tech program of china. h.l. was also supported by knowledge innovation program of the chinese academy of sciences.
for three different kinds of binding sites, mosfom displays an excellent ability to differentiate active compounds through energy and shape complementarity. the multi-objective optimization method was successfully applied in virtual screening with two different scoring functions that can yield reasonable binding poses and can furthermore, be ranked with the potentially compromised conformations of each compound, abandoning those conformations that can not satisfy overall objective functions. different scoring functions focus on diverse aspects of ligand binding, and no single scoring can satisfy the peculiarities of each target system. results suggest that mosfom can effectively enhance enrichment and performance compared with a single score. in this paper, two multi-objective optimization methods, called mosfom, were developed for virtual screening, which simultaneously consider both the energy score and the contact score. therefore, the idea of a consensus score strategy was put forward.
 development of a fast and accurate scoring function in virtual screening remains a hot issue in current computer-aided drug research. different scoring functions focus on diverse aspects of ligand binding, and no single scoring can satisfy the peculiarities of each target system. therefore, the idea of a consensus score strategy was put forward. integrating several scoring functions, consensus score re-assesses the docked conformations using a primary scoring function. however, it is not really robust and efficient from the perspective of optimization. furthermore, to date, the majority of available methods are still based on single objective optimization design.  in this paper, two multi-objective optimization methods, called mosfom, were developed for virtual screening, which simultaneously consider both the energy score and the contact score. results suggest that mosfom can effectively enhance enrichment and performance compared with a single score. for three different kinds of binding sites, mosfom displays an excellent ability to differentiate active compounds through energy and shape complementarity. efmoga performed particularly well in the top 2% of database for all three cases, whereas moea_nrg and moea_cnt performed better than the corresponding individual scoring functions if the appropriate type of binding site was selected.  the multi-objective optimization method was successfully applied in virtual screening with two different scoring functions that can yield reasonable binding poses and can furthermore, be ranked with the potentially compromised conformations of each compound, abandoning those conformations that can not satisfy overall objective functions. 