Quantitative clustering and matching of conformers using average electron densities

1 Introduction

Interactions among molecules are highly dependent on the conformation of the structures (Lundqvist et al., 2004; Matsumori et al., 2005; Roach et al., 2005; Atkins et al., 2001; DeLorbe et al., 2009). In drug design, for instance, conformations play a crucial role in the “key & lock” complementarity (Basciu et al., 2019; Verdonk et al., 2008; Tripathi and Bankaitis, 2018; Chiarparin et al., 2019). This complementarity impacts, among others, the binding affinity, the biological activity, and the potency of a drug (Kastritis and Bonvin, 2013; Cera, 2020; Kenny, 2019; Cournia et al., 2017). This strong dependence of the “key & lock” complementarity on the conformational changes of molecules urges the need for a tool that can classify conformers based on how they can interact with a given receptor (or any other molecule).

ESP maps are 3D maps calculated from charge densities, taking into account the electron density and the nuclear charge distribution of a molecule. They visually illustrate the electron-rich and electron-poor areas by the displaying positive and negative lobes (Francisco et al., 2021; Bolcato et al., 2022; Arabi, 2020; Matta et al., 1868), which highlight the electrophilic and nucleophilic regions in a molecule (Bolcato et al., 2022). Hence, ESP maps are commonly used in the literature to guide the assessment of the complementarity between a molecule and its receptor (Matta, 2014). Molecules with similar ESP maps exhibit comparable interactions or complementarities with a given receptor. Although it is possible to visually group conformers based on the similarity in their ESP maps, this task is prone to subjectivity, particularly when dealing with a large number of conformers. It is rather difficult, and even challenging, to avoid ambiguities in grouping ESP maps by comparing their shapes, especially that these maps are highly dependent on the minor variations in the chosen isodensity values of the maps (see Fig. 1).

Close

Fig. 1

Electrostatic potential maps of ibuprofen at isodensity values of 0.25, 0.50, 0.75, 1.00, and 1.25 atomic units (au). Red and blue lobes represent regions of negative and positive electrostatic potentials, respectively.

Export to PPT

In the literature, there are approaches that cluster conformers based on structural parameters (Gupta and Payne, 2001; Takasu et al., 2000; Kessel et al., 2001), or based on energies (Sakhaee et al., 2016; De et al., 2017). However, these clustering approaches are not helpful for the task of assessing the “key & lock” complementarity with a given receptor. This is because (as will be shown in the results and discussion) similarities in structures do not necessarily correlate with similarities in ESP maps. Similarly, conformers with similar energies do not necessarily have comparable ESP maps. Thus, it is desired to have a tool that meticulously and robustly cluster conformers in groups that share similar ESP maps, without ambiguities. Such a tool can facilitate the computational part of drug discovery. For example, given that all conformers within one group share similar ESP maps, indicating comparable “key & lock” complementarity with a given receptor, it is possible to assess the docking of only one conformer per group, as a representative, instead of having to test all conformers. This could save time and computational resources.

This paper demonstrates that the AED tool is capable of accurately clustering conformers into groups that share comparable ESP maps. The clustering is based on similarities in the AED values of a particular group of interest within a molecule, e.g. a bioisostere within a drug molecule. To evaluate the AED of the atoms constituting the group of interest within a molecule, it is important to be able to partition the molecule into atomic basins. This can be done using a partitioning scheme such as QTAIM (Bader, 1990; Matta and Boyd, 2007). According to the QTAIM theory, zero-flux surfaces partition the molecular space into atomic basins. In other words, an atom in a molecule is bound by a surface that has zero-flux in the gradient vector field of the electron density. For outer atoms where one of the atomic basin sides extends to infinity, this side is delimited by an outer isodensity envelope. Once the atomic basins are defined, atomic properties including, among others, atomic electron populations, charges, volumes, and AED values can be evaluated. The AED is the ratio of electron population to volume (Arabi, 2020; Matta et al., 1868; Matta, 2014; Arabi, 2017; Arabi and Matta, 2016) of an atom, a group of atoms within a molecule, or an entire molecule. It is given by $${\mathit{AED}}_{\mathit{group}}=\frac{{\sum }_i^{\mathit{al}l\mathit{atoms}\mathit{in}\mathit{the}\mathit{group}}N(\mathrm{\Omega })}{{\sum }_i^{\mathit{all}\mathit{atoms}\mathit{in}\mathit{the}\mathit{group}}V(\mathrm{\Omega })},$$ where N(Ω) is the electron population per atom Ω, and V(Ω) is the atomic volume. For example, the AED of the carboxylic acid group (–COOH) in ibuprofen is equal to the sum of the atomic electron populations of the carbon atom, two oxygen atoms, and the hydrogen atom in the –COOH group, divided by the sum of the atomic volumes of the same atoms within this group.

In this paper, conformers of ibuprofen will be quantitatively grouped according to the AED of its carboxylic acid group. In another example, conformers of a tetrazole analogue of ibuprofen (see Fig. 2) will be grouped according to the AED values of its tetrazole group. The success of the clustering will be confirmed by inspecting the similarities in the ESP maps of the conformers within each group. This is the first time that the AED tool is used in the context of clustering conformers of a given molecule into groups that share similar ESP maps.

Close

Fig. 2

2D structures of ibuprofen (a) and its tetrazole analogue (b), along with their respective molecular properties obtained from the OpenEye Software.

Export to PPT

Another aim of this study is to assess the capability of this AED tool to match conformers of different molecules. This will be tested via the concept of bioisosterism in drug design. Bioisosteric substitutions in drug molecules are performed to fine-tune the pharmacokinetic and pharmacodynamic properties of lead compounds. Non-classical bioisosteric moieties have significant variations in their populations, charges, shapes, polarities, number of H-bond donors and acceptors, etc. (Arabi, 2020; Matta et al., 1868; Arabi, 2017; Arabi and Matta, 2016; Arabi, 2021; Matta and Arabi, 2011), yet the lead compound and its substituted analogue maintain similar biological activities (Meanwell, 2011; Meanwell, 2013; Brown, 2012; Langmuir, 1919). This similarity is believed to stem from similarities in their ESP maps and, thus, their “key & lock” complementarity with their receptor. A typical example is the substitution of a carboxylic acid group with a tetrazole group (Brown, 2012; Ballatore et al., 2013; Birch et al., 2009; Zou et al., 2020). Ibuprofen is a non-steroidal anti-inflammatory (NSAID) drug that contains a carboxylic acid group. Ibuprofen was recently suggested as treatment for type II diabetes, fatty liver disease, and atherosclerosis through inhibiting the human adipocyte fatty-acid binding protein (FABP4) (González and Fisher, 2015). To optimize the properties of ibuprofen in this drug repurposing application, its carboxylic acid bioisostere could be replaced with a tetrazole group. The AED tool has been previously explored in the context of bioisosterism. It was shown that capped bioisosteres have similar AED values (Arabi, 2020; Matta et al., 1868; Arabi, 2017; Arabi and Matta, 2016; Arabi, 2021; Matta and Arabi, 2011; Osman and Arabi, 2022). In particular, the carboxylic acid and tetrazole bioisosteres, when capped with a methyl group, have AED values that differ by no more than 0.2 % (Matta et al., 1868). However, in a later study (Osman and Arabi, 2024), and based on the findings of the first part of this study, it was observed that conformational changes in larger molecules (e.g. ibuprofen) lead to variations in the AED values of the bioisosteric moiety. These variations among conformers are in the range of 2 %, i.e. one order of magnitude greater than the 0.2 % AED difference observed between the carboxylic acid and tetrazole bioisosteres. Thus, in this study, the AED tool will also be tested, for the first time, for its capability to match conformers of different molecules that are meant to exhibit similar “key & lock” complementarity within a given receptor.

2 Computational methods

Conformers of the (S)-ibuprofen, (S)-2-[4-(2-methylpropyl)phenyl]propanoic acid, and its tetrazole analogue, 5-[1-[4-(2-methylpropyl)phenyl]ethyl]-2H-tetrazole, are considered in this study (see Fig. 2).

To estimate the deformation of the receptor structure by the presence of the (S)-ibuprofen ligand, the FABP4/(S)-ibuprofen complex (PDB ID: 3P6H) was superimposed with the apo FABP4 protein (PDB ID: 3RZY) (González and Fisher, 2015). Using the Pymol package (Schrödinger, 2015), the RMSD was found to be 0.51 Å, which is about 4 times smaller than the acceptable RMSD values reported in references (Ramírez and Caballero, 2018; Hao et al., 2015) and references therein. The (S)-ibuprofen coordinates were extracted from the PDB structure. The Omega toolkit (Hawkins et al., 2010) in OpenEye was used to generate all conformers of (S)-ibuprofen (43 conformers) and its tetrazole analogue (40 conformers). The Make Receptor tool was used to prepare the receptor (https://www.eyesopen.com) in a 37.09 Å × 40.20 Å × 50.01 Å box that includes the solvating water molecules. All conformers of both ligands were docked into their receptor using the FRED tool (McGann, 2012) within the OpenEye Software. VIDA 4.4.0.4. Santa Fe, NM. (https://www.eyesopen.com) was used to view the electrostatic potential maps.

Single point calculations were performed on all conformers, using the Gaussian16 package (Frisch et al., 2016). The simulations were completed using the same method used in previous studies to evaluate AED values (Arabi, 2020; Matta et al., 1868; Arabi, 2017; Arabi and Matta, 2016; Osman and Arabi, 2022), that is the B3LYP functional in density functional theory (Becke, 1993), with the 6–311++G(d,p) basis set (Pople et al., 1989), ultrafine grids, and tight SFC criteria.

Using the AIMALL software (Keith, 2014), the QTAIM analysis was performed on the wavefunction files of the quantum simulations. The outer isodensity envelope of the atomic basins is taken as the 0.001 atomic units (au) envelope. This isodensity covers more than 99 % of the electron population of the atom in a molecule, and the volumes at this isodensity are comparable to the experimental van der Waals volumes in the gas phase (Bader, 1990; Bader, 2001). The maximum Lagrangian value was 4.4 × 10⁻⁴ au.

The AED values were evaluated for the carboxylic acid and tetrazole groups in ibuprofen and its tetrazole analogue, respectively. These groups were selected based on two reasons. In the first and second parts of the study, the clustering of the conformers of a given molecule was based on the AED value of its most electronegative group, i.e. the carboxylic acid and tetrazole groups in ibuprofen and its analogue, respectively. In the third part of this study, the matching between conformers of different molecules is tested based on the bioisosteric connection between their carboxylic acid and tetrazole groups.

SPSS (version 28) was used to perform K-means clustering of the AED values, where K was set to 11 for both ibuprofen and its analogue.

3 Results and discussion

3.1

3.1 Part I: Clustering of the conformers of ibuprofen

Fig. 3A shows the K-means grouping of the ibuprofen conformers based on the AED values of the carboxylic acid functional group. Fig. 3B depicts the similarity in the ESP maps of the clustered conformers within each group. The grouping is accurate in all cases, except for conformer 36. This conformer is, in fact, an outlier according to the K-means analysis. In K-means clustering, the farther a point is from the cluster center (which is the average of all the data points in this particular cluster), the less likely it is for this point to belong to that cluster. The Euclidean distance between conformer 36 and the cluster center of G6 (i.e. the group it belongs to) is 5.0 × 10⁻⁵ au. This is, by far, the largest among all cases, it is almost 66.7 % greater than the second largest distance, 3.0 × 10⁻⁵ au (for conformer 41 in G3).

Close

Fig. 3

(A) AED of the carboxylic acid group in 43 conformers of ibuprofen. The conformers are clustered in 11 groups. Below each ESP map is the corresponding conformer number, as depicted in Part A. (B) ESP maps of the clustered conformers. Red and blue lobes represent regions of negative and positive electrostatic potentials, respectively.

Export to PPT

It is worth noting that the full range of AED values spans a maximum difference of only 2.26 %. This reflects the robust power of this AED tool in meticulously clustering the conformers despite the small differences in the AED values of their bioisosteric moieties. Fig. 3B clearly shows that all conformers within the same group share almost identical ESP maps, with similar shapes and sizes of the positive and negative lobes. Some of the groups, e.g. G8 and G9, may seem to have similar ESP maps, but upon close inspection, it is visually noticeable that the separations between the upper and lower red lobes in the ESP maps of G8 versus G9 are different. In another example, G9 and G10 may seem to have same ESP maps, but the distribution of their red lobes is different. The ESP maps in G6 and G11 have a big red lobe that is almost planar in G6, but rather out of plane in G11. G2 and G6 also appear to have similar ESP maps, but the angle of the ring plane (relative to the carboxylic acid moiety) is different between the two groups. Moreover, in G2, the big red lobe extends slightly wider, while in G6, it is relatively thicker. This non-obvious meticulous grouping also highlights the power of the AED tool at depicting the minor differences in the ESP maps. The close similarities in the AED values and ESP maps per group is unique and cannot be accurately detected by visual inspection of the ESP maps. Even the 2D shapes (or energies) of the molecules do not correlate with the similarities in the ESP maps per group. For example, conformers 10 and 20 of ibuprofen (shown in Fig. 3B) have different 3D structures, yet they have similar ESP maps and they belong to the same group, G2.

3.2

3.2 Part II: Clustering of the conformers of the tetrazole analogue of ibuprofen

Fig. 4A shows the K-means clustering based on the AED values of the tetrazole group in 40 conformers of the tetrazole analogue of ibuprofen. The precision in this clustering is remarkable given that the AED values of the bioisosteric moiety in all conformers span a small range from 0.0709 to 0.0725 au (i.e. a maximum difference of 3.50 %). The ESP maps of ibuprofen (Fig. 3B) have bigger red negative lobes and smaller blue positive lobes than those of the tetrazole analogue (Fig. 4B). It is obvious from Fig. 4B that the clustering is precise. Two exceptions can be noted: (i) conformer 17 in G8 and (ii) G5 could be split into two subgroups. Similar to the case of ibuprofen, some groups, e.g. G7 and G9 here, may seem to have similar ESP maps. However, a careful comparison shows that the orientation of the ring plane (relative to the tetrazole group) and its two associated small red lobes on either side of the ring is different in the two groups. Additional examples of similar groups could be mentioned, but they would have differences in size, orientation, and/or arrangement of the lobes.

Close

Fig. 4

(A) AED of the tetrazole group in 40 conformers of the tetrazole analogue. The conformers are clustered in 11 groups. Below each ESP map is the corresponding conformer number, as depicted in Part A. (B) ESP maps of the clustered conformers. Red and blue lobes represent regions of negative and positive electrostatic potentials, respectively.

Export to PPT

3.3

3.3 Part III: Drug design application: matching groups of conformers of ibuprofen and its tetrazole analogue

All conformers of ibuprofen and its tetrazole analogue were docked into their receptor. The top-ranked conformers of ibuprofen and its tetrazole analogue have similar Chemgauss4 docking scores of −10.53 and −10.14, respectively. This difference in the decimal places is rather irrelevant given the approximate nature of the docking scores. Therefore, as expected biologically from the definition of bioisosterism, and as noted chemically from the similar docking scores, the bioisostere-containing ibuprofen and its tetrazole analogue share similar “key & lock” complementarity with their receptor. This similarity is further emphasized by the large percent similarity of 98.77 % in the average AEDs of the bioisosteric groups (see Table 1). Table 1 shows the percent similarities in the AED values of the carboxylic acid and the tetrazole bioisosteric moieties among various conformers of ibuprofen and its tetrazole analogue, respectively. It is obvious that the percent AED similarities listed in this table are very high, exceeding 95 % in all comparisons reported. This highlights the thoroughness of the AED tool in depicting AED similarities among the carboxylic acid and tetrazole bioisosteric moieties, despite the conformational changes in their respective molecules. In addition, the small AED difference of 1.23 % in the average AED values of the bioisosteric groups is contrasted with the large AED differences observed between non-bioisosteric groups, i.e. groups that are not supposed to share similar “key & lock” complementarity with a given receptor. For example, the average AED of the carboxylic acid bioisostere in ibuprofen (0.0725 ± 5.95 × 10⁻⁴ au) is 26.48 % greater than the average AED of the rest of the ibuprofen molecule (0.0533 ± 2.28 × 10⁻⁴ au), noting that these two parts of the molecules are not bioisosteres of each other. Similarly, the average AED of the tetrazole bioisosteric moiety in the tetrazole analogue (0.0716 ± 4.51 × 10⁻⁴ au) is 25.70 % greater than the AED of the rest of the molecule (0.0532 ± 1.73 × 10⁻⁴ au), again noting that these two parts of the molecule are not bioisosteres of each other.

Table 1 Percent similarities in the AED values of the bioisosteric moieties in selected conformers or groups of conformers of ibuprofen and its tetrazole analogue.

Conformer of ibuprofen	AED (in au) of the carboxylic acid moiety	Conformer of the tetrazole analogue of ibuprofen	AED (in au) of the tetrazole moiety	Percent similarity in the AED values (%) of the bioisosteric moieties
The conformer with the minimum AED value	0.0715	The conformer with the minimum AED value	0.0709	99.19
The conformer with the maximum AED value	0.0740	The conformer with the maximum AED value	0.0725	97.97
The average of all conformers	0.0725	The Average of all conformers	0.0716	98.77
The conformer with the maximum AED value	0.0740	The conformer with the minimum AED value	0.0709	95.87
The top-ranked conformer in docking, (conformer 1)	0.0721	The top-ranked conformer in docking, (conformer 1)	0.0714	99.04
The top-ranked conformer in docking, (conformer 1)	0.0721	The third-ranked conformer in docking, (conformer 3)	0.0721	99.98
The average of the conformers in G2 from Fig. 3	0.0720	The average of the conformers in G8 from Fig. 4	0.0720	99.98
The average of the conformers in G3 from Fig. 3	0.0721	The average of the conformers in G9 from Fig. 4	0.0721	99.97

Fig. 5A shows that the AED values for the bioisosteres in the conformers of G2 of ibuprofen and G8 of its tetrazole analogue are similar, and so are the respective ESP maps in these two groups (Fig. 5B). Similar observations are made with respect to G3 (associated with ibuprofen) and G9 (associated with the tetrazole analogue). The docking results align with those obtained from the AED tool. The top-ranked ibuprofen (conformer 1, Chemgauss4 docking score of −10.1 units) belongs to G3 and the third-ranked tetrazole analogue (conformer 3, Chemgauss4 docking score of −10.0 units) belongs to G9; noting that these two groups are detected with the AED tool to share similarities (see Fig. 5). The AEDs of the bioisosteric moieties in these two conformers are 99.98 % similar, which is more than the percent AED similarity observed for the two top-ranked conformers (99.04 %). This 99.98% similarity is ca. 4 % greater than the smallest percent AED similarity; i.e. 95.87 % between conformer 13 of ibuprofen and conformer 2 of its analogue (see Table 1). It is worth noting that the first- and third-ranked tetrazole analogues (conformers 1 and 3, respectively) share similar ESP maps, especially when comparing the blue lobe and the two red lobes around the tetrazole moiety; the difference is only in the angle of the two red lobes on both sides of the benzene ring of the molecule (see Fig. 4). It is obvious, from Fig. 3B and Fig. 4B, that the ESP of the top-ranked conformers of ibuprofen and its tetrazole analogue share similarities in the positive and negative lobes and in their rough positions. However, the similarity is ambiguous enough not to deterministically decide on whether or not the two molecules are likely to have similar interactions with their receptor. On the other hand, the AED difference between the bioisosteric moieties of the top-ranked ibuprofen and its tetrazole analogue is 0.96 %. This difference of less than unity is deterministic and remarkable, especially given the noticeable and numerous differences in the molecular properties of ibuprofen and its tetrazole analogue as summarized in Fig. 2).

Close

Fig. 5

(A) Matching of the AED values between the groups of conformers of ibuprofen and its tetrazole analogue. (B) ESP maps of the matched groups of conformers. Below each ESP map is the conformer number as depicted in Fig. 3 for ibuprofen and Fig4 for the tetrazole analogue of ibuprofen. Red and blue lobes represent regions of negative and positive electrostatic potentials, respectively.

Export to PPT

Overall, despite the differences in the two bioisosteric moieties (see Fig. 2), and the numerous possible conformer pairs among the 43 ibuprofen conformers and 40 tetrazole analogue conformers, the AED tool accurately matches conformers of different molecules, with similarities up to 98.77 % (see Table 1).

4 Conclusions

In conclusion, the AED is a unique tool used to quantitatively and meticulously cluster conformers of a molecule into groups that share similar ESP maps and, therefore, share similarities in their “key & lock” complementarity with a given receptor. This AED tool is reliable and also has a remarkable accuracy exceeding 96 %. This AED tool is powerful not only in clustering conformers of a single molecule, but also in matching conformers of different molecules in a manner that highlights shared similarities in their ESP maps. The ability of the AED tool to identify molecular structures with analogous ESP maps help expedite drug discovery and optimization processes.

Declaration

Two patent applications related to this work were filed to the United States Patent and Trademark Office. One is granted (Method of classifying conformers, patent number: 11862295, Inventor: Alya A. Arabi, Assignee: United Arab Emirates University, Date of patent: January 2, 2024), and the other one is still under review.

Credit authorship contribution statement

Alya A. Arabi: Conceptualization, Funding acquisition, Data curation, Writing – original draft, Writing – review & editing, Visualization, Investigation, Validation, Formal analysis, Methodology, Resources, Project administration, Software.