The shape similarity of two compounds is calculated as follows: (1) VABg=∑i∈A,j∈BVijg=∑i∈A,j∈BpipjKij(παi+αj)32, where pi and pj are set to 2.7, αi and αj obtain the van der Waals value for each atom and (2) Kij=exp(−αiαjRij2αi+αj), where Rij is the distance between atoms i and j.

Notice that the accuracy obtained from (1) depends on the number of atoms in the two compared molecules, i.e. the higher this number, the longer the value of VAB as an absolute value. To be able to measure the level of similarity between compounds, regardless of the number of atoms that they are composed of and the descriptor used, the Tanimoto Similarity (Jaccard, 1901) value is computed as follows: (3) Tcs=VABVAA+VBB−VAB, where VAB is the A molecule overlaid onto B molecule. VAA and VBB is the overlap of the molecules A and B, respectively. (3) has a value in the range [0,1] , where 0 means there is no overlapping and 1 means the shape of both molecules is the same.

The electrostatic similarities are obtained by numerical solution of the Poisson equation (Böttcher et al., 1974), viz: (4) ∇{ϵ(r)∇ϕ(r)}=−ρmol(r), where ϕ(r) is the electrostatic potential, ϵ(r) is the dielectric constant, and ρmol(r) is the molecular charge distribution. Electrostatic similarity between two compounds is compared by determining EAB : (5) EAB=∫ϕA(r)ϕB(r)ΘA(r)ΘB(r)dr≈h3∑ijkϕijkAϕijkBΘijkAΘijkB, where Θ is a masking function to ensure potentials interior to the compound are not considered part of the comparison. The integral appearing in (5) is a volume integral, computed using a grid-spacing parameter, h.

Again the accuracy obtained by (5) depends on the number of atoms in the compared molecules. As such, similarly to what was done previously, the Tanimoto Similarity (Jaccard, 1901) value has been computed as follows: (6) TcE=EABEAA+EBB−EAB, where EAB is the A molecule overlaid onto B molecule. EAA and EBB is the overlap of the molecules A and B, respectively. In this case, (6) has a value in the range [−0.33,1] , where −0.33 means the charges of both compounds have the same value but opposite loads, 0 means there is no overlapping, and 1 means the charges of both molecules are the same.

All the experiments in this work have been executed using a Bullx R424-E3, which consists of 2 Intel Xeon E5 2650v2 (16 cores), 128 GB of RAM memory and 1 TB HDD (http://hpca.ual.es/en/infraestructure) along with the cluster Eagle https://wiki.man.poznan.pl/hpc/index.php?title=Eagle.

In this work, a database provided by The Food and Drug Administration has been used (FDA). The Food and Drug Administration is a federal agency of the United States Department of Health and Human Services responsible for protecting and promoting public health by controlling, among other things, prescription and over-the-counter pharmaceutical drugs (medications). This agency provides a data set containing 1751 compounds, which represents approved medicines that can be safely used on humans in the USA. This database is useful since in the high similarity cases it would directly contribute to drug re-purposing. This is of relevant utility given the clear trend regarding re-purposing drugs observed over the last 5 years (Dakshanamurthy et al., 2012; Kumar and Zhang, 2018; Yuan et al., 2017).

The version of the database used in this work was obtained from DrugBank v5.0.1 (Wishart et al., 2018) and necessary mol2 files for the VS calculations were set up by using AmberTools (Case et al., 2017) by removing salts and neutralizing their protonation state, computing partial charges by MMFF94 force field, adding hydrogen atoms and minimizing energies (default parameters) (Halgren, 1995).

A comprehensive computational analysis may cover a representative sample of the database. The compounds included in the FDA database have different attributes, one of the most relevant for the study at hand being the number of atoms. In this work, a selection of 50 compounds has been made in the following way: the compounds in the database have been sorted by the number of atoms, including hydrogen, and then divided into 24 intervals (see Fig. 2). From each sector, at least one compound was chosen at random and proportional to the number of compounds in the sector.

Finally, these comparisons between compounds have been run using OptiPharm_ES with the following input parameter configuration: N=200000 function evaluations, M=5 starting poses, tmax=5 iterations and Rtmax=1 as the smallest possible radius.

As previously mentioned, the LBVS-Shape bases its predictions on a pre-selection of the first best compounds in terms of superimposition score (N). In this subsection, a study has been conducted to know how the value of N affects the final results from the point of view of electrostatic similarity. In particular, the LBVS-Shape has been performed on the selected 50 queries and for five different values of N, i.e. N has been set to 175, 438, 876, 1313 and 1751 compounds. It means that for each query, we have selected either 10% , 25% , 50% , 75% or 100% of the ranked compounds during the pre-selection phase.

Figure 3 illustrates a toy example of the main steps of the LBVS-Shape method for the Query DB01213 and N=1751 , i.e. the total number of compounds in the FDA set. Initially, the Query is compared to each compound TargetS from the database to obtain their optimum position and corresponding shape similarity value TcS . As previously mentioned, this stage is carried out by using ROCS. Afterwards, compounds are sorted ( RkS ) in decreasing order by TcS . The N best compounds are selected and evaluated to measure the corresponding electrostatic similarity value TcEEval . Notice that the evaluation of the electrostatic similarity considers the pose obtained with the shape similarity optimization. The compound with the highest TcEEval , called BestComp throughout this paper, is selected as the best prediction. Finally, as an additional and unconsidered stage in the LBVS-Shape method, we have computed the optimized superposition between the BestComp and the Query by using OptiPharm_ES. The corresponding TcE value is then provided.

To get an overview of the results, average values of the BestComp found for the 50 queries and each value for N have been computed, and shown in Table 1. In particular, the average position Av(RkS) in the sorted list where the BestComp were located have been computed, together with the following: their mean number of atoms Av(NS) , their average shape similarity value Av(TcS) , their corresponding electrostatic similarity value Av(TcEEval) when they are evaluated, and finally, their mean electrostatic similarity when they are optimized Av(TcE) .

As it can be seen, the predictions seem to improve in term of electrostatic similarity as the number N of selected molecules in the sorted list increases (see columns Av(TcEEval) and Av(TcE) ). In accordance with these results, the posterior comparison between LBVS-Shape and LBVS-Electrostatic methods has been carried out by setting N=1751 .

To analyse the performance of both methods, we have conducted a study in which the selected 50 molecular queries are processed with reference to the FDA database. Notice that comparing a query with itself always reaches the maximum similarity value, both for electrostatic potential as well as for shape. Subsequently, these results were removed when ranking the compounds. In other words, the compounds given as a result are not the most similar ones, but the second compounds in the ranked list. Additionally, as previously mentioned, the traditional method has been carried out considering the total number of compounds in the database N=1751 , so as to increase the probability of finding better predictions.

To illustrate how we generate the later summarizing tables, a sample of the results obtained by both methods when comparing a query to the molecules in the dataset is studied. In particular, the instance Query=DB01213 is analysed. Notice that this is the example used to illustrate the stages of the LBVS-Shape method in Fig. 3. After that, the same instance is considered to exemplify the performance of the LBVS-Electrostatic method (see Fig. 4). Notice that this Query has been selected because it is small and it helps to see the main ideas of the paper very easily by using figures. However, the conclusions inferred from the associated results can be extrapolated to any other Query . As can be observed, the LBVS-Electrostatic technique solves an optimization problem to determine the electrostatic similarity, TcE , between the pharmaceutical Query and every TargetE in the database. Afterwards, the list of compounds is sorted by the TcE value and the one located in first position, RkE=1 , is selected as the best prediction. Finally, to complete the study, optimization is carried out to calculate the shape similarity TcS between the chosen compound and the Query .

For the sake of clarity and comparison, the results shown in Figs. 3 and 4 are summarized in Table 2. The meaning of the columns as well as the particular values in the tables, are the ones previously explained and shown in each figure. The last row corresponds to the values associated with the best predictions. As can be observed, each method obtains a different compound as a top solution. LBVS-Shape provides the DB00184 molecule with a TcS=0.621 and a TcEEval=0.500 . At the same time, LBVS-Electrostatic proposes the DB03255 compound as being the most similar to the query with TcE=0.810 and TcSEval=0.880 . As such, the LBVS-Electrostatic method has not only obtained a more similar compound in terms of electrostatic potential, but also in shape. In Fig. 5, the final position for each case is shown.

Once the specific case of DB01213 has been explained in detail, the results of the 50 queries have been summarized in Table 3. Columns RkEEval and RkE have been removed in this table because their values are always 1. The last row summarizes the average of the results.

As evidenced, LBVS-Electrostatic obtains on average TcE=0.738 , which is higher than that given by LBVS-Shape, TcEEval=0.497 . Similar conclusions can be inferred when comparing the TcE average values for both methods. Additionally, when the results are analysed individually, we can see that LBVS-Electrostatic provides solutions with higher TcE values than those achieved by LBVS-Shape. In fact, in 48 out of 50 cases, LBVS-Electrostatic obtains a different compound than that reached by LBVS-Shape.

Regarding shape similarity, it is possible to infer that, on average, the methods are equivalent in terms of accuracy of the predictions, i.e. LBVS-Shape obtains an average value of Tcs=0.554 while LBVS-Electrostatic reaches a mean value of Tcs=0.505 . Furthermore, analysing the obtained results individually, we can see that in 2 out of 50 cases, LBVS-Electrostatic offers better or equivalent predictions than that achieved by LBVS-Shape in terms of shape (see columns Tcs in LBVS-Shape and TCsEval in LBVS-Electrostatic). It means that cases exist where two compounds can be very similar in terms of electrostatic potential, although they can be very different in terms of three-dimensional shape. It means that those solutions could not be obtained by using the methodology followed by the traditional LBVS-Shape method, since it only focuses on the compounds with the highest similarity in shape.

Making a somewhat more detailed approach for compounds smaller than 50 atoms, which means the first 23 query compounds in the table, there are 5 cases where the difference is less than 0.05 (DB00529, DB00173, DB00331, DB00915 and DB01352) and in another 3 cases the difference is 0.1 (DB01043, DB07615 and DB01268). Considering the values of these 7 cases in which the shape LBVS-Electrostatic is smaller than that of LBVS-Shape, the average difference is 0.048, while the mean gain in electrostatic similarity for those 7 compounds is 0.271. In large compounds, which includes 27 queries, there are only two cases with similar characteristics, which are compounds DB09236 with a difference of 0.07 and DB06699 with a difference of 0.013, both of them for shape similarity. In view of these results, the LBVS-Electrostatic method seems to be justified when proposing new solutions for small compounds.

However, not all the improvements are related to electrostatic fields. The optimization of electrostatic potential using OptiPharm_ES might allow a better solution to be found in terms of shape too. Compounds DB01119 and DB1213 in Table 3 are some outstanding examples. For example, in the case of Query=DB01119 , the best compound found by LBVS-Shape is DB00828 with TcS=0.655 and TcEEval=0.519 . Moreover, LBVS-Electrostatic’s best compound is DB00173. It has a better TcE , i.e. 0.829, but also the position of those compounds after the electrostatic optimization is improved, TcSEval=0.779 .

The ZAP Toolkit has been widely used in the literature to calculate the electrostatic similarity score for two compounds (Boström et al., 2013; Tresadern et al., 2009; Chu and Gochin, 2013; Kim et al., 2015; Kossmann et al., 2016; Woodring et al., 2017; Maccari et al., 2011; Kim et al., 2016; López-Ramos and Perruccio, 2010; Hevener et al., 2012; Kaoud et al., 2012; Tiikkainen et al., 2009; Massarotti et al., 2014; Oyarzabal et al., 2009; Haque and Pande).

In this subsection we would like to remark that the ZAP Toolkit can return an erroneous value, which was discovered when using OptiPharm_ES. During the optimization procedure, OptiPharm_ES can progressively separate two input compounds aimed to escape from local optima and explore the searching space in depth. In fact, it is possible to analyse cases where no overlap exists between the input molecules. During the analysis of the results, we discovered that cases exist where the ZAP Toolkit can overflow, mainly when situations such as the previously mentioned happen. See Fig. 6 to see a particular example. Herein, compound DB01365 remains fixed on the left while compound DB00459 occupies three positions (red, blue and pink). The red compound obtains an electrostatic similarity value of 1. The light blue compound is displaced half a unit to the left, i.e. closer to the reference compound and its similarity value is 0.38. The pink compound is shifted 0.5 units to the right, that is, away from the reference compound. Its similarity value is 0. Calculations can be made using the ZAP Python script available at https://docs.eyesopen.com/toolkits/python/zaptk/thewayofzap.html in the Electrostatic Similarity section.

This problem has been solved in OptiPharm_ES by considering the poses with the previously mentioned problem unfeasible. It means that they are no longer considered during the optimization process.