This paper aims to maximize both the shape T c S and the electrostatic similarity T c E of two molecules (a query A and a target B), simultaneously.

The shape similarity score is calculated using a weighted Gaussian model (Grant and Pickup, 1995; Grant et al., 1996; Yan et al., 2013), which is widely used in other works (Yan et al., 2013; Cui et al., 2015; Lo et al., 2016) and can be expressed mathematically as follows: (1) V A B = ∑ i ∈ A , j ∈ B w i w j v i j , where v i j is the sum of the intersections of the atoms and w i and w j are weights associated with the atoms i and j, respectively. From this function, the Tanimoto coefficient (Real and Vargas, 1996) is computed to normalize the score: (2) T c S = V A B V A A + V B B − V A B , where V A A and V B B are the overlaps of the molecules A and B with themselves, respectively. This function returns a value in the range [ 0 , 1 ] where 0 means that there is no overlap and 1 means that there is a complete overlap between both molecules.

The electrostatic similarity (Böttcher and Bordewijk, 1978) between A and B is obtained by calculating the (3): (3) E A B = ∫ ϕ A ( r ) ϕ B ( r ) Θ A ( r ) Θ B ( r ) dr ≈ h 3 ∑ i j k ϕ i j k A ϕ i j k B Θ i j k A Θ i j k B .

In a similar fashion to (1), the Tanimoto metric is used to normalize the value. For (4), the range value is [ − 0.33 , 1 ] where 0 means there is no electrostatic overlap, 1 means the overlap is complete and − 0.33 means the overlap is complete but the loads of the molecules are opposite. (4) T c E = E A B E A A + E B B − E A B .

For optimization purposes, we look for the target’s pose that maximizes both the shape and electrostatic score values. To do so, we consider that the query remains in the same position during the whole optimization process. A pose p represents the rotation and translation of the target from its original position, and it is defined as a quaternion of the form p = { θ , c 1 , c 2 , Δ }, where θ is the angle that a molecule is rotated with respect to the axis created by the points c 1 = ( c 1 x , c 1 y , c 1 z ) and c 2 = ( c 2 x , c 2 y , c 2 z ). Δ = ( Δ x , Δ y , Δ z ) is a displacement vector.

The multi-objective optimization problem can be then defined as: (5) T c ( p ) = { max T c S ( p ) , max T c E ( p ) } , s.t. p ∈ S .

The need to tackle this problem from a multi-objective perspective can be seen in Fig. 1. It shows an example where the query molecule DB01155 and the target DB01208 are considered as input. Figure 1(a) shows the pose obtained when the shape similarity is maximized. As can be seen, a high value of T c S = 0.832 is obtained. However, when we evaluate such a pose with the electrostatic similarity, we realize the obtained value is quite small ( T c E = 0.556). In Fig. 1(b) we show the opposite case, i.e. we optimize the electrostatic similarity and evaluate the shape similarity with the obtained pose. As can be seen, we obtain high values of T c E = 0.892 but a negligible value for shape similarity T c S = 0.318. All this clearly means that there exist instances where both objectives are in conflict and then, multi-objective techniques are needed.

When multiple objectives are present, the concept of an optimal solution, as in the single-objective problems, does not exist (Pardalos et al., 2017). Then, before addressing the optimization problem, it is necessary to explain what solving a multi-objective problem means. In a single-objective context it is easy to determine when a pose p is better than another p ′ . It clearly happens when T c i ( p ) > T c i ( p ′ )T{c_{i}}({\boldsymbol{p}^{\prime }})$]]>, being i ∈ { S , E }. In a multi-objective framework, the quality of two poses is determined by the concept of dominance. Its definitions are as follows:

Figure 2 illustrates the concept of dominance. Notice how the objective space is represented, i.e. we have mapped a set of solutions p to the objective space T c ( p ) = ( T c S , T c E ).

Therefore, solving a multi-objective problem like the one in (5) means finding the whole non-dominated subset formed by all the efficient poses, whose corresponding T c = ( T c S , T c E ) represents the optimal Pareto front (Deb et al., 2016). However, obtaining an exact description of the efficient set (or Pareto front) is impossible for the problem at hand since those sets are continuous and include an infinite number of points. Furthermore, the computing cost may be high, which is an important point in the current context. As such, in this work, we will focus on providing a finite set of points comprising a Pareto Front Approximation (PFA) as a solution of (5) (see Fig. 3).

This subsection describes the methodology used in the literature that is characterized by the use of single-objective software for the selection of candidate compounds. As will be seen, regardless of the case, not only is a significant effort required to determine which molecules to select for evaluation, but the selection is highly dependent on the available information and expert knowledge, which can rule out promising results.

The most widely used methodology in the literature consists of facing the query molecule to the whole database by optimizing the shape similarity. This procedure, denoted as Mono-Shape throughout this paper, is illustrated through an example in Fig. 4(a), where the Query molecule is DB09074, and the database consists of the rigid compounds from the well-known FDA. Initially, the Query is compared to each compound Target S from the database to obtain their optimum position and corresponding shape similarity value T c S . Afterward, a post-processing procedure starts, whose complexity and required effort depends on the experience of the decision-maker. In this respect, there can be many different workflows. As part of this work, three of them will be explained.

The first workflow ( W F 1 Mono-Shape ) is trivial and consists of the decision-maker only being interested in similarity of shape. In such a case, the compounds are sorted ( R k S ) in descending order by T c S , and the molecule with the greatest T c S is proposed as the best prediction. In our example, this molecule would be DB11799 with T c S = 0.69.

However, if the decision-maker’s preferences include molecules with both high shape and electrostatic similarities, which we named the second workflow ( W F 2 Mono-Shape ), the post-processing is rather time-consuming. This workflow is widely used in the literature and for more details on the specific procedures, several papers can be consulted (Chu and Gochin, 2013; Kim et al., 2015; Kossmann et al., 2016; Woodring et al., 2017). After sorting the molecules by the T c s value, it continues by selecting and evaluating the N best compounds to measure the corresponding electrostatic similarity value T c E Eval . Notice that the evaluation of the electrostatic similarity considers the pose obtained with the shape similarity optimization. The decision-maker usually determines the number N of analysed molecules and usually is not greater than 10% of the database sizes (Hevener et al., 2012; Kaoud et al., 2012; Kossmann et al., 2016). Even so, the number of compounds to be analysed might be considerably high, meaning the selection of promising solutions might demand a lot of experience, knowledge, and time from the expert’s side. In our example, for N = 100 and after a visual inspection of the aligned poses, the expert could select the molecule DB13801 with a T c S = 0.48 and T c E Eval = 0.32 grounded on several criteria (additional chemical substructural similarities and potential known pharmacophores for the chosen drug target, among others).

Finally, another possible way to proceed might be to compute the average value between T c S and T c E Eval and propose the one with the greatest mean value ( T c meanSE ) as the best prediction, which is called the third workflow ( W F 3 Mono-Shape ). It is intended to alleviate a problem that frequently occurs in W F 2 Mono-Shape . As part of this, compounds can have a very high value in shape similarity but very low in electrostatic similarity.

Let us now itemize the obtained predictions when we optimize the electrostatic similarity instead of the shape similarity, denoted as Mono-Elec. We can now consider the expert has a single-objective tool available to maximize the electrostatic similarity of two molecules. As previously, to illustrate the procedure, the same query and database are selected (see Fig. 4(b)). Similarly, the query molecule faced the whole database by optimizing the electrostatic similarity of each target, and a list of candidate molecules sorted according to their T c E values is obtained. From this point on, the experts can make different decisions based on their knowledge, available information, or preferences. They might consider that only a descriptor is of interest and hence select the molecule with the greatest T c E , ( W F 1 Mono-Elec ). In our example, this molecule would be DB00956 with a T c E = 0.48. They might also require molecules with high similarities in both descriptors. As previously considered, the experts might select the top N molecules and evaluate the obtained pose with the shape similarity descriptor, to ultimately select the one with the greatest evaluated value, ( W F 2 Mono-Elec ). In our example, for N = 100 the expert would propose DB14840 with T c E = 0.29 and T c S Eval = 0.47. Either that, or they might compute the average value ( W F 3 Mono-Elec ) and propose the molecule DB11656 as the best prediction with T c meanES = 0.39.

To summarize both methodologies, Mono-Shape and Mono-Elec, it has been shown that the decision-maker needs a great deal of knowledge and experience to be able to carry out many decisions ( W F 1 , W F 2 , W F 3 ), solve different single-objective optimization problems and accomplish different analyses in search of good predictions. The key point is that all of them may lead to different candidate molecules that could be of interest to the experts. However, the larger the number of candidate molecules obtained, the bigger the number of optimization problems to solve, and the greater the post-processing effort required by the decision-makers, meaning it becomes unfeasible for a human expert to apply their experience and chemical knowledge to such large scale contexts.

This section describes the methodology proposed in this work. For this purpose, it has been divided into two parts. Firstly, we will detail the tool, called MultiPharm, which has been developed to support this methodology. Secondly, MultiPharm-DT will be explained. It is the compound selection strategy based on a decision tool that implements the expert’s preferences.

The execution of MultiPharm over a database will obtain a Pareto front for each pair (query, target). It could be too much information for the expert, even considerably more than the one managed with the single-objective perspective. In this work, we have also implemented a decision tool, called MultiPharm-DT, to reduce the quantity of information given to the expert. Although one could consider many different selection criteria, we opted to implement a mechanism that obtains the Pareto front of the paretos as part of MultiPharm-DT.

Figure 7 helps to illustrate this mechanism. As can be seen, MultiPharm-DT merges all the particular Pareto fronts. At this point, the experts can compare them, see the similarities of the compounds from the database with the query molecule, and extract information in different ways. Apart from the union of the paretos, MultiPharm-DT applies a decision mechanism and provides a filtered solution. The expert is able to select the most appropriate compound for their preferences at a glance. That means they will select the compound DB00956 if they are looking for a molecule with high electrostatic similarity. If their preferences are inclined towards shape similarity, they will choose the compound DB11799. However, if they are interested in proposing compounds with a balance between the two descriptors, they will choose from the compounds DB01012, DB01224, DB01259 or DB11656. Finally notice that for each compound, the multi-objective methodology offers a set of possible poses, as mentioned in subsection The Multi-Objective Optimization Problem, which could also be of great interest to the decision-maker.

The database used in this work was obtained from DrugBank v5.0.1 (Wishart et al., 2018) and mol2 files necessary for the VS calculations were set up by using AmberTools (Case et al., 2017). This took the form of removing salts and neutralizing their protonation state, computing partial charges by MMFF94 force field, adding hydrogen atoms and minimizing energies (default parameters) (Halgren, 1995) (see Fig. 8). From this database, we selected those compounds that were rigid and validated by the FDA, a federal agency of the United States Department of Health and Human Services. This subset was classified by its number of atoms in groups. A compound (query) was selected randomly from each group, as shown in Fig. 9. Consequently, 32 molecules were selected for analysis against the rest of the database. And in order to compare software and methodologies on an equal footing, hydrogen atoms were considered for shape and electrostatic similarity.

All the experiments in this work have been executed using 18x Bullx R424-E3, which consists of 2 Intel Xeon E5 2650 (16 cores), 64 GB of RAM memory and 128 GB SSD (http://hpca.ual.es/en/infraestructure).

To compare the two methodologies, we decided to use OptiPharm in its two versions (Puertas-Martín et al., 2019; Puertas-Martín et al., 2020) and MultiPharm. OptiPharm is a recent piece of software designed to work with the LBVS problem which implements an evolutionary global optimization algorithm that can solve any optimization problem related to the similarity of two compounds, a query and target molecules. It implements procedures to adjust the latter to the query molecule, remaining unchanged throughout the optimization process. OptiPharm has proved to be very competitive when maximizing the shape similarity, as well as the electrostatic. For a more comprehensive explanation of the software and the obtained results, reading the original works (Puertas-Martín et al., 2019; Puertas-Martín et al., 2020) is recommended. It is also available in the free-to-use server BRUSELAS (Banegas-Luna et al., 2019).

OptiPharm and MultiPharm can be parameterized in order to adapt them to different problems and the preferences of the user. Focusing on the quality of the solutions, the configuration established for OptiPharm is the same as published in other works. This will be called OptiPharm Robust (OpR) and its input parameters are: N = 200,000 function evaluations, M = 5 starting poses, t max = 5 iterations and l max = 1 as the smallest possible radius. Conversely, MultiPharm is set up with the default parameters (Zhang and Li, 2007; Li and Zhang, 2008). But in order to compare both pieces of software on equal terms, the number of MultiPharm evaluations is set to the same value as the OptiPharm configuration, i.e. the number of evaluations is 200,000. In addition, the number of solutions that will form a Pareto front is set at 300.

In this subsection, an exhaustive analysis of the results from Table 2 is performed, showing the values of T c S and T c E that could not be displayed in the table due to the large number of poses offered by MultiPharm-DT. Thus, 3 groups can be distinguished. In the first group there are 9 queries characterized by a high similarity value in both descriptors. In addition, MultiPharm-DT finds the same compounds as W F 1 Mono-Shape and W F 1 Mono-Elec . Two examples are the DB00536 and DB01242 queries. The most similar target for the first one has similarity values of T c S = 0.967 and T c E = 0.800 while regarding the target for the second query, the values are T c S = 0.987 and T c E = 0.919 (Fig. 12).

The second group of compounds consists of 14 queries. It is characterized by MultiPharm finding more compounds in addition to those from W F 1 Mono-Shape and W F 1 Mono-Elec . Figure 13 shows an example of 6 different queries. These queries are of great interest and are where the potential of MultiPharm-DT can truly be seen. Normally, the compounds found by W F 1 Mono-Shape and W F 1 Mono-Elec are very good in terms of the descriptor they optimize but not the other. However, with MultiPharm, in addition to getting these two compounds, it is possible to find a set where the balance of both descriptors is much better and by giving a slightly lower value to one descriptor, the similarity of the other improves considerably.

Finally, the third group consists of the remaining 9 queries, where MultiPharm finds several molecules, but among them, only one belongs either to W F 1 Mono-Shape or W F 1 Mono-Elec . Four of these cases are represented in Fig. 14 and in fact, in Fig. 14(c) it can be seen that among the MultiPharm compounds, the DB13801 molecule, and not the DB00427 molecule, is found by W F 1 Mono-Elec .

In view of all the above results, it can be seen that this new methodology with MultiPharm-DT is able to widen the knowledge about the compounds to be selected since different poses with similarity values in both descriptors are known for each compound. Indeed, it improves the quality of the decisions and allows good compounds to be selected without the knowledge or action of a decision-maker, which greatly facilitates the appropriate choice of compounds.

Finally, as a general comment, it is important to mention that MultiPharm-DT provides more valuable predictions than the single-objective workflows, in terms of quantity and quality. Interestingly, this is in a single run without the expert’s participation. Furthermore, without considering the expert’s effort, obtaining the solutions for the single-objective methodology will cost 200,000 evaluations per workflow (a total of 400,000 evaluations) while MultiPharm-DT only requires 200,000. It is clear to see then that the advantages of using MultiPharm-DT, in terms of computing time, are sizeable.