Article Cite This: J. Chem. Inf. Model. 2018, 58, 165−181
pubs.acs.org/jcim
xMaPAn Interpretable Alignment-Free Four-Dimensional Quantitative Structure−Activity Relationship Technique Based on Molecular Surface Properties and Conformer Ensembles Jan Dreher,‡ Josef Scheiber,# Nikolaus Stiefl,† and Knut Baumann* Institute of Medicinal and Pharmaceutical Chemistry, University of Technology Braunschweig, Beethovenstrasse 55, D 38106 Braunschweig, Germany S Supporting Information *
ABSTRACT: A novel alignment-free molecular descriptor called xMaP (flexible MaP descriptor) is introduced. The descriptor is the advancement of the previously published translationally and rotationally invariant threedimensional (3D) descriptor MaP (mapping property distributions onto the molecular surface) to the fourth dimension (4D). In addition to MaP, xMaP is independent of the chosen starting conformation of the encoded molecules and is therefore entirely alignment-free. This is achieved by using ensembles of conformers, which are generated by conformational searches. This step of the procedure is similar to Hopfinger’s 4D quantitative structure−activity relationship (QSAR). A five-step procedure is used to compute the xMaP descriptor. First, a conformational search for each molecule is carried out. Next, for each of the conformers an approximation to the molecular surface with equally distributed surface points is computed. Third, molecular properties are projected onto this surface. Fourth, areas of identical properties are clustered to so-called patches. Fifth, the spatial distribution of the patches is converted into an alignment-free descriptor that is based on the entire conformer ensemble. The resulting descriptor can be interpreted by superimposing the most important descriptor variables and the molecules of the data set. The most important descriptor variables are identified with chemometric regression tools. The novel descriptor was applied to several benchmark data sets and was compared to other descriptors and QSAR techniques comprising a binary fingerprint, a topological pharmacophore descriptor (Cats2D), and the field-based 3D-QSAR technique GRID/PLS which is alignment-dependent. The use of conformer ensembles renders xMaP very robust. It turns out that xMaP performs very well on (almost) all data sets and that the statistical results are comparable to GRID/PLS. In addition to that, xMaP can also be used to efficiently visualize the derived quantitative structure−activity relationships.
■
INTRODUCTION Structure−activity correlation techniques need to represent molecules numerically. During the past years and decades a considerable number of molecular descriptors have been developed. Programs such as Dragon1,2 or VCCLab3 compute up to 5000 descriptors per molecule. These descriptors can be distinguished by the way they represent molecules. Onedimensional descriptors (1D) include bulk parameters as well as physicochemical properties (e.g., log P, molecular volume). Two-dimensional (2D) techniques use information on the molecular graph. Three-dimensional (3D) descriptors are based on molecular geometry (i.e., the Cartesian coordinates of the molecule’s atoms) and can roughly be divided into two different groups: grid-based techniques such as CoMFA,4 CoMSiA,5 GRID,6 or Continuous Molecular Fields,7 which transform grids into continuous functions, and distance-based techniques such as MaP,8,9 the Fuzzy Pharmacophores,10 and GRIND.11 Both types use a different frame of reference. In the grid-based techniques changes of the molecule’s position in space cause changes in the descriptor (external frame of reference). The position of the molecule in space is irrelevant if distances between molecular features are used to characterize it © 2017 American Chemical Society
(internal frame of reference) since distances between two or more molecular features of a single conformer (i.e., a rigid object) do not change when the object is moved in space (translational and rotational invariance). If an external frame of reference is used, the molecules need to be aligned in space to derive meaningful descriptors. This latter step is obviously not necessary for distance-based techniques and thus avoids potential bias caused by the chosen alignment rule. The downside of this approach is a more difficult visualization and interpretation of the resulting model. However, for several newer distance-based approaches this problem has been solved.8−11 In 3D quantitative structure−activity relationship (QSAR), only a single conformer per molecule is analyzed. Hence, an appropriate conformer for each molecule in the data set needs to be chosen. Since most biologically active molecules are more or less flexible, the latter selection step represents another potential source of bias. Moreover, with a single conformer 3DQSAR methods cannot represent molecular flexibility. This has Received: July 13, 2017 Published: November 27, 2017 165
DOI: 10.1021/acs.jcim.7b00419 J. Chem. Inf. Model. 2018, 58, 165−181
Article
Journal of Chemical Information and Modeling
The article is organized as follows. In the next section the method is described in detail and compared further to other approaches. In the results part of the study, the novel molecular structure descriptor is benchmarked against alignment independent 2D-QSAR methods and alignment dependent 3D-QSAR methods using 12 well-established data sets, eight of which originate from Sutherland et al.40 Two-dimensional descriptors are included as simple approaches to check whether there is a real benefit in terms of predictivity of the more sophisticated methods. In addition to the benchmark study, xMaP is studied in detail for acetylcholinesterase inhibitors.
early been realized by Hopfinger and co-workers who introduced the first 4D-QSAR technique which encodes molecular geometry plus molecular flexibility.12 A lot of publications show the power of this particular 4D-QSAR technique.13−17 Briefly, Hopfinger’s 4D-QSAR uses a gridbased method (external frame of reference) and requires an alignment of the molecules prior to analysis. It uses large, homogeneous conformer ensembles generated by molecular dynamics simulations. As a grid-based method, interpretation of the 4D-QSAR models is straightforward. Similar to Hopfinger’s 4D-QSAR, the newly introduced xMaP approach is also based on conformer ensembles to reflect molecular flexibility. As opposed to Hopfinger’s 4D-QSAR, it uses an internal frame of reference and thus needs no alignment. The conformers are generated by conformational searching and cover a larger conformational space. Overall, xMaP is quite different in its nature resulting in different strengths (no alignment). Dobler and Vedani extended receptor-surface models18 to encode molecular flexibility.19,20 This technique has also been extended to the fifth and sixth dimensions representing target flexibility and different solvation states.21,22 It matured over the years and is broadly applicable. However, as a receptor-surface model, the technique is also alignment dependent. Another alignment dependent multiconformational method related to receptor-surface models was published recently.23−25 Similar to Hopfinger’s 4D-QSAR, homogeneous conformer ensembles are employed here. Self-organizing maps (SOM) have also been used successfully for 4D-QSAR.26−29 Here the coordinates and the partial charges of the conformational ensemble of each molecule are used to construct a 2D-SOM which is then used as input for PLS regression with variable selection. Again, this method requires the alignment of the molecules since the SOM construction uses an external frame of reference. LQTA-QSAR in essence extended CoMFA to the fourth dimension and also is alignment dependent.30,31 It uses molecular dynamics simulations to generate the conformational ensemble and a grid-based approach with different probes to generate CoMFAlike descriptors. An enormous number of descriptors (several ten thousands) is generated so that variable filtering and variable selection becomes very important.32 An alignment independent approach that uses multiple conformers was developed by Kuz’min and colleagues.33−35 Here, a set of different simplexes (predefined tetratomic structural fragments with fixed composition, chirality and symmetry) is used to encode the conformers. This approach to 4D-QSAR does not take long-range distances into account because simplexes consist of one sphere of neighboring atoms. Moreover, it results in a huge number of variables. Another 4DQSAR method extends the electron-conformational (EC) method,36 which is based on quantum chemistry, to conformational ensembles37 and also is alignment independent. Since the size of the so-called electron-conformational matrices of conjunction (ECMC) depends on the number of atoms, the method first performs a pharmacophore elucidation to identify identically composed fragments of constant size (so-called electron conformational submatrices of activity) to enable QSAR modeling. The method is thus restricted to data sets of closely related analogs with a common scaffold. A recent approach also uses molecular dynamics successfully to characterize the conformational space.38 It uses alignmentindependent WHIM descriptors to encode the molecules.39 While these descriptors holistically characterize each conformer, they cannot easily be back-projected onto the molecule.
■
MATERIALS AND METHODS xMaPOverview. The initial step of the newly developed technique is the computation of conformer ensembles for each molecule in the data set by a conformational search procedure. If the 3D structures within one ensemble form more than one cluster, only the largest one is kept (this will be referred to as harmonization of a conformer ensemble). That means, that only a single binding mode is modeled which is chosen to be the energetically most likely one. Afterward a discretized molecular surface is computed for each of the remaining conformers. Each of the surface points is assigned a hydrophobicity property (strongly hydrophilic, weakly hydrophilic, strongly hydrophobic, weakly hydrophobic). In addition, surface points can be assigned a H-bonding property (H-bond donor, H-bond acceptor) if applicable. Next, the large number of surface points is reduced to surface patches by clustering areas of identical properties. These surface patches are represented by their centers of mass and surface size. They are then used to generate potential two-point-pharmacophores. These are characterized by a particular patch property combination and the distance between the centers of mass of the two patches. The resolution on the distance axis is 1 Å. Put differently, each vector entry stores the occurrence of a particular two-point-pharmacophore for a range spanning 1 Å (e.g., 7.0 ≤ x < 8.0 Å). This information is initially stored in a single vector for each conformer. Once all conformers are encoded, the mean value for the occurrence of each potential two-point-pharmacophore is computed across all conformers. The resulting mean vector is used to represent a single molecule. The computation of the descriptor is depicted in Figure 1 and explained in more detail in the following. (i) Determination of the Conformer Ensemble. In xMaP, the fourth dimension is described by the conformational states a molecule can adopt. Initially, molecular dynamics (MD) simulations were used to sample the conformational space in accordance with previous publications by Hopfinger and coworkers.41 These authors use an upper energy limit of 2 kcal/ mol over the minimum structure and calculate 50 000 conformers per molecule. In this procedure only the energetically most favorable conformational space is sampled. By default, xMaP uses conformational searches and also samples the more heterogeneous less favorable conformational space (up to 25 kcal/mol) which was found to be important in order to reproduce the bioactive conformation.42 In addition to the different coverage of the conformational space, MD simulations are computationally more expensive than is standard conformational searching. For example, when modeling the PGF2α-data set shown below, the same statistical quality with MD simulations and conformational searching was obtained at the cost of a 10-fold increase in computation time (data not shown). 166
DOI: 10.1021/acs.jcim.7b00419 J. Chem. Inf. Model. 2018, 58, 165−181
Article
Journal of Chemical Information and Modeling
Where NA and NB is the number of conformers in cluster CA and CB respectively. The term d(ai, bj) is the average Euclidean distance between conformer ai ∈ CA and bj ∈ CB d(ai , bj) =
1 Natoms
Natoms
∑ k=1
(xai − xbj )2 + (ya − yb )2 + (zai − z bj )2 k
k
ik
jk
k
k
(2)
Where Natoms is the number of heavy atoms and x, y, z represent the respective coordinates. At each stage of the hierarchical clustering, the clusters CA and CB with minimum D(CA, CB) are merged. The merging stops if the minimum distance is equal to or greater than D(CA, CB) ≥ 0.45 Å (default). The conformers of the largest cluster are kept which are thought to represent the conformational space of the energetically most favorable (potential) binding mode. Elimination of outlying confomers from the ensemble might seem inappropriate for a 4D technique that intends to address molecular flexibility. However, harmonization primarily eases the visual interpretation of the derived QSAR equation. The statistical quality of the results is only marginally affected (see Results). As was pointed out during the reviewing process, encoding a conformer ensemble using a SOM (self-organizing map)47 as in SOM4D-QSAR26,48 may lead to a similar result as the clustering described above since SOMs reinforce the respective output neurons for similar conformers. (iv) Calculation of the Molecular Surface and Mapping of the Molecular Properties. For each remaining conformer in the data set, a discretized molecular surface consisting of surface points is computed. The calculation follows MaP and was previously described.8 Briefly, the surface points are calculated by a modified GEPOL algorithm49 with a grid spacing of 0.8 Å. The comparatively coarse grid spacing is insensitive to the position of the molecule in the grid owing to the formation of surface patches (see below). Subsequently, the molecular lipophilic potential (MLP, i.e. hydrophilicity/hydrophobicity)18,50 is computed for each surface point. Depending on the actual MLP a category is assigned according to the scheme shown in Figure 2.
Figure 1. Calculation of the xMaP descriptor: (i) Compute conformer ensemble. (ii) Superimpose conformers to the lowest energy conformer. (iii) Harmonize conformer ensemble, i.e. remove outlying conformers. (iv) Determine discretized molecular surface. (v) Cluster areas of identical properties to surface patches. (vi) Encode patch combinations as potential two-point-pharmacophores and store their occurrence.
In xMaP conformational searching is performed by employing Omega43 from OpenEye Inc. This conformational search algorithm was extensively evaluated.42,44 The HQS (highquality screening) setting42 was used as default setting for xMaP. It uses an energy window of 25 kcal/mol above the energy minimum of a particular structure and a maximum number of 500 conformers. Conformers with an RMSD value smaller than 0.4 Å to conformers already in the ensemble were discarded. Identical parameters were used for all data sets. (ii) Superimposition and (iii) Harmonization of the Conformers. For each conformer ensemble the structure with the lowest (predicted) energy is determined and serves as template. Next, each remaining conformer is superimposed onto the template by employing the Kabsch algorithm.45,46 It is important to stress that this is not an alignment step in the sense of a 3D-QSAR analysis because conformers of a single molecule are fit to a reference of the same molecule. This results in a conformer ensemble with minimum overall RMSD. In order to eliminate outlying conformers, which are far away from the ensemble mean, average-linkage clustering is performed. The distance between two clusters is computed as the average distance between conformers of the first cluster to the conformers of the second cluster. Mathematically, this distance D(CA, CB) between two clusters CA and CB (the linkage function) is computed as N
D(CA , CB) =
Figure 2. Classification of surface points according to the molecular lipophilic potential (MLP).
Apart from the MLP, surface points associated with atoms that show H-bonding potential are additionally assigned either the H-bond donor or the H-bond acceptor property. This is in contrast to MaP where each surface point could show only a single property which is disadvantageous here since a noncontiguous surface with respect to the MLP would result. Positive and negative ionized groups are encoded here through an increase in hydrophilicity of the surface and changes in the H-bonding pattern. As was pointed out by one reviewer, a more direct encoding of charges might ease interpretation of the results and will be considered in future versions of the descriptor. Encoding charges implicitly at present is mainly done to keep the number of variables appreciably smaller.
N
A B 1 ∑ ∑ d(ai , bj); NA ·NB i = 1 j = 1
ai ∈ CA , bj ∈ CB (1) 167
DOI: 10.1021/acs.jcim.7b00419 J. Chem. Inf. Model. 2018, 58, 165−181
Article
Journal of Chemical Information and Modeling (v) Generation of Surface Patches. In the next step, contiguous surface areas exhibiting identical properties are merged to form so-called patches. Surface patches for potential H-bond donors and H-bond acceptors always consist of the surface area of the particular atom which shows the H-bonding potential (i.e., for a hydroxyl group, the surface area of the hydrogen forms the H-bond donor patch). The MLP categories are grouped to patches by applying a combination of single and average linkage clustering. First, single linkage clustering on the 3D coordinates of a particular MLP category is performed using the SLINK implementation51 to identify contiguous surface points (each member of this subset has at least one neighbor within s · 3 , where s represents the grid spacing; default: s = 0.8 Å). These contiguous surface points form a patch with ns_slink points. If ns_slink exceeds a cutoff, the patch is either split into several patches if it is large enough or it is reduced to a smaller one which has the effect of shifting the center of mass closer to the functional group that gives rise to this particular patch (see below). In both cases, average linkage clustering is performed using the ns_slink surface points. Here, the distance between two clusters, that define potential patches, is defined as the average distance between all pairs (a, b) of surface points, where a is a surface point from the first cluster and b is a surface point from the second cluster (cf. eq 1). The distance d(a, b) between two surface points a ∈ CA and b ∈ CB is the Euclidean distance (cf. eq 2). NA and NB is the number of surface points in clusters A and B, respectively. Surface points are grouped together until D(CA, CB) reaches 6 Å. This cutoff resulted from visual inspection of many surface patches and represents a value where a patch usually covers an entire functional group. If the size of the thus generated patches does not exceed a cutoff of minimal surface points, they are deleted. Average linkage clustering ensures that only densely populated surface properties (concerning a particular property) are considered. Sparsely populated border regions are deleted. This smoothing is mainly done to find the correct center of mass for the bulk of the data. Once this step is completed, the molecular surface is divided into regions of different properties. With this information each conformer can be represented by the distribution of its patches (cf. Figure 1). Using the patches instead of the single surface points to characterize the molecule leads to a massive reduction in data, distills the essential information to form potential two-point-pharmacophores, and eases the obligatory interpretation of the resulting models. To characterize the importance of a patch, the surface area A is computed and stored for each patch as A = ns ·s 2
distribution is stored in histograms. For the six pharmacophoric properties used here (H-bond acceptor (A), H-bond donor (D), weakly lipophilic (Lw), strongly lipophilic (Ls), weakly hydrophilic (Hw), strongly hydrophilic (Hs)), 21 different types of P2PPs result which are named according to the two properties involved (e.g., LsA). Their occurrence is stored in 21 different histograms which are chained in a prespecified order to form a vector that characterizes the encoded conformer. The default bin width in the histogram is 1 Å. The distance range starts at 0 Å and ends at distance dMax, which depends on the data set. Plain histograms count the number of objects that fall into its categories. Here, each entry in the histogram is weighted by the product of the underlying patch sizes W (where W = Ai·Aj for patches i and j). Hence, these weighted histograms resemble Moreau and Broto’s autocorrelation and cross-correlation vectors52 with patch centers instead of atomic coordinates and patch sizes instead of atomic properties. In order to avoid discretization errors due to the arbitrary bin borders and thus large fluctuations of the descriptor for very similar conformers, not only the matching bin of the histogram is incremented but also the neighboring bin (so-called fuzzy increment53). The increment for each of the two bins depends on the difference of the actual patch−patch distance to the center of the matching bin. The increment of the matching bin is computed as follows ⎛ | bin − d | ⎞ ⎟ inc = W · ⎜1 − ⎝ binW ⎠
(4)
where d denotes the distance, bin is the center of the matching bin, and binW the bin width. If bin − d < 0, then the neighboring bin with a larger distance is incremented by W − inc, else the bin with the next smaller distance is incremented. After all P2PPs have been recorded, the resulting vector of chained histograms describes a specific conformer by the spatial distribution of its surface patches weighted by patch size. The total number of variables is m = 21dMax, i.e. the different types of P2PPs times the number of distance bins. The value of dMax is set to the distance of the most distant patch pair rounded to nearest integer greater than this particular distance plus one bin. This last bin is added to collect the fuzzy increment of the penultimate bin. The encoding is repeated for all conformers of the molecule under scrutiny. The result is a matrix of size nc × m (nc = number of conformers, m = number of variables) where each row describes a single conformer and each column describes the weighted occurrence of a defined P2PP over all conformers. For further computations this matrix is transformed to the ensemble average descriptor (i.e., take the column mean). Conformational flexibility is captured as follows: For patches in flexible parts of the molecule the weight for the associated P2PPs will not accumulate in a single distance bin but will be spread out over several bins. This distribution may not be confused with fuzzy increments but is due to different distances between the patches in different conformers. Patches in rigid parts of the molecules will give rise to P2PPs that increment predominantly the same distance bin(s). That way the resulting average vector characterizes both, the distribution of surface properties, and the flexibility of the molecule. This is illustrated in Figure 3 where three different descriptor vectors for a single molecule are shown. The dashed and the dotted lines characterize two different conformations. These descriptor vectors are very different, since their underlying conformations are extremes within the raw (i.e., not yet harmonized)
(3)
where ns is the number of surface points forming the patch and s the grid spacing. If a patch was smoothed by means of average linkage clustering, the importance of the original patch is kept by using ns_slink instead of ns to compute the surface area. It should be recalled here that patch smoothing is primarily done to find the proper center of mass and not to change the size of the patch. (vi) Calculation of the Descriptor. The actual molecular descriptor characterizes the distribution of the molecule’s potential two-point-pharmacophores (P2PP). A P2PP is characterized by two patches and the Euclidean distance separating the centers of mass of these patches according to the general scheme (patch property 1)−distance−(patch property 2). All P2PPs of a particular conformer are computed and their 168
DOI: 10.1021/acs.jcim.7b00419 J. Chem. Inf. Model. 2018, 58, 165−181
Article
Journal of Chemical Information and Modeling
Table 1. Summary of the Applied Validation Protocol 1
2
3 4 5
Figure 3. Illustrative example. Descriptor vectors of two conformations (dashed and dotted) and the ensemble average vector (bold).
6
conformer ensemble. The bold line represents the average vector over the conformer ensemble. Here, flexible parts of the molecule are encoded by smooth, broad peaks because the underlying patch−patch distances differ a lot from conformer to conformer. Model Generation and Validation. Not all of the P2PPs will be related to the biological activity under scrutiny. Hence, an important statistical tool for the novel descriptor is variable selection. By and large, model building and model validation follows previously established protocols54−56 with an important modification which assesses test set predictivity more stringent in a two-layered validation scheme using either bagging or leave-multiple-out cross-validation in the outer loop to generate test set partitions.57 A brief overview of the steps used in model generation and validation can be found in Table 1. For variable selection, the reverse-elimination-method tabusearch (REM-TS) with an early stopping rule is employed.58 The search is stopped when changes in the objective function from one iteration to the next iteration were only small (
