Probing the origins of human acetylcholinesterase inhibition via QSAR modeling and molecular docking

Data set

A data set of inhibitors against human AChE (Target ID CHEMBL220) were compiled from the ChEMBL 20 database (Gaulton et al., 2012) that is comprised of a total number of 9,242 bioactivity data points from 5,049 compounds. SMILES notations of the compounds were curated with ChemAxon’s Standardizer (ChemAxon Kft., 2015) using the same parameter settings as described in our previous study (Simeon et al., 2016). The initial data set was assembled from several bioactivity measurement units including (in order of decreasing data size) IC₅₀, K_i, % activity, % inhibition, MIC, EC₅₀, etc. IC₅₀ was selected for further investigation as they constituted the largest subset with 4,910 compounds. A closer look revealed that 1,301 compounds had no reported IC₅₀ values or had lesser/greater than signs which were subjected to removal thereby leaving 3,609 compounds. Only data points having nM as the bioactivity unit were selected for further study, which produced 3,596 compounds. Furthermore, redundant compounds having different bioactivity values were kept if the standard deviation of IC₅₀ was less than 2 and this resulted in 2,571 compounds. Moreover, some compounds were found to have no SMILES notation associated with it and were thus removed. A final data set comprising of 2,570 compounds was obtained.

Description of inhibitors

AChE inhibitors were encoded by a vector of fingerprint descriptors accounting for its molecular constituents. Prior to calculating descriptors, salts were removed and tautomers were standardized using the built-in function of the PaDEL-Descriptor software (Yap, 2011).

Although fingerprint descriptors are able to capture the feature space of chemical compounds, their ability to be used as descriptors for bioactivity modeling can vary. In fact, performance differences existing amongst the different fingerprint type has been the subject of several investigations into its utilization for bioactivity modeling. Riniker & Landrum (2013) benchmarked and assessed the performance of predictive models constructed from 2D fingerprint descriptors obtained from RDKit.

In this study, the suitability of 12 different fingerprint descriptors for predicting the bioactivity of AChE inhibitors was investigated. Table 1 summarizes the employed fingerprints along with their corresponding size, description and reference.

Additionally, the four molecular descriptors that are used to define the Lipinski’s rule-of-five comprising of molecular weight (MW), logarithm of the octanol/water partition coefficient (ALogP), number of hydrogen bond donor (nHBDon) and number of hydrogen bond acceptor (nHBAcc) were also computed by the PaDEL-Descriptor software.

Feature selection

Collinearity is a condition where descriptor pairs are known to have intercorrelation, which not only add complexity to the model but could potentially give rise to bias. To remedy this, the cor function from the R package stats was used to find the pairwise correlation among descriptors, and descriptors in a pair with a Pearson’s correlation coefficient greater than the threshold of 0.7 was filtered out using the findCorrelation function from the R package caret to obtain a smaller subset of descriptors (Kuhn, 2008).

Data splitting

To avoid the possibility of bias that may arise from a single data split when building predictive models (Puzyn et al., 2011), predictive models were constructed from 100 independent data splits and the mean and standard deviation values of statistical parameters were reported. The data set was split into internal and external sets in which the former comprises 80% whereas the latter constitutes 20% of the initial data set. The sample function from the R base package was used to split the data.

Multivariate analysis

Supervised learning is to learn a model from labeled training data which can be used to make prediction about unseen or future data (James et al., 2013). This study constructs regression models, which affords the prediction of the continuous response variable (i.e., pIC₅₀) as a function of predictors (i.e., fingerprint descriptors).

Random forest (RF) is an ensemble classifier that is composed of several decision trees (Breiman, 2001). Briefly, the main idea behind RF is that instead of building a deep decision tree with an ever-growing number of nodes, which may be at risk for overfitting and overtraining of the data, rather multiple trees are generated as to minimize the variance instead of maximizing the accuracy. As such, the results will be more noisier when compared to a well-trained decision tree, yet these results are usually reliable and robust. The ranger function from the R package ranger, which is a fast implementation of the RF algorithm that was used for constructing the models (Wright & Ziegler, 2015).

Validation of QSAR models

Model validation is an important process, which should be performed to ensure that a fitted model can accurately predict responses for future or unknown subjects. Two statistical parameters were used to evaluate the performance of the QSAR models consisting of Pearson’s correlation coefficient (r) and root mean squared error (RMSE). The r value is a commonly used metric to represent the degree of relationship between two variables of interest. It can range from −1 to +1 in which negative values are indicative of negative correlation between two variables and vice versa. RMSE is a commonly used parameter to assess the relative error of the predictive model. The predictive performance of the QSAR models was verified by 10-fold cross-validation, external validation and Y-scrambling test.

The 10-fold cross-validation technique does not used the entire data set to build predictive model. Instead, it splits the data into training and testing data set by allowing model that is built with training data set us allow to assess the performance of the model on the testing data set. By performing repeats of the 10-fold validation, the average accuracies can be used to truly assess the performance of the predictive model.

Y-scrambling test was used to ensure the robustness of the predictive model not only to rule out the possibility of chance correlations but also to assess the statistical significance of R² and Q², ensuring the generalizability of QSAR model. The true Y-dependent variable (i.e., pIC₅₀) was randomly scrambled and the statistical assessment parameters are recalculated. Performance of the Y-scrambling test can be deduced from the regression line of the plot of R² versus Q². Intercept values for R² and Q² as denoted by iR² and iQ², respectively, were calculated. Negative iQ² is indicative of an acceptable QSAR model and that there is no chance correlation from the real model (Eriksson et al., 2003). Furthermore, $r_m^2$ and $\mathrm{\Delta }{r}_m^2$ metrics as introduced by Roy et al. (2013) were used to verify the robustness of the proposed QSAR model in which an acceptable QSAR model should give $r_m^2>0.5$ and $\mathrm{\Delta }{r}_m^2<0.2$.

Applicability domain analysis

The applicability domain (AD) estimates the likelihood of reliable prediction for compounds on the basis of how similar they are to compounds used to build the model. Thus, compounds falling outside the AD may lead to unreliable predictions. The most common approach for determining AD is described by Gramatica (2007) and Tropsha, Gramatica & Gombar (2003), which is to compute the leverage values for each compound. The leverage value allows one to identify whether new compounds will lie within or outside the domain. Leverage values for all compounds are calculated via adjustment of X to give the hat matrix H: (1)${\displaystyle H=X{\left(X^{T}X\right)}^{-1}{X}^T}$

where X is a two-dimensional matrix comprising of n compounds and m descriptors while X^T is the transpose of X. Meanwhile, the leverage value of the ith compound (h_i) is the ith diagonal element of H: (2)${\displaystyle h_i=x_i^T{\left(X^{T}X\right)}^{-1}{x}_i}$

where x_i is the descriptor row-vector of the ith compound. The warning leverage h^∗ is calculated by: (3)${\displaystyle h^{\ast }=3\left(p+1\right)∕n.}$

Practically, the leverage value along with the William’s plot is often used to assess the AD of QSAR models. The William’s plot is constructed by depicting the standardized residuals versus the leverage value for each compound’s h_i. If the ith compound has h_i > h^∗ then it means that the ith compound exerts a great influence on the QSAR model and may be excluded from the AD. In spite of this, it does not appear to be an outlier because its standardized residual may be small.

Molecular docking

The co-crystal structure of human AChE with donepezil (PDB ID: 4EY7) was retrieved from the Protein Data Bank and initially prepared by removing alternative side chains and water molecules. The protein was prepared via the rebuilding of bonds and the addition of missing hydrogen atoms. Subsequently, the protein was cleaned by merging the atomic charge and removing lone pair atoms, non-polar hydrogen atoms and non-standard amino acid residues. Grid box was set up to provide coverage of the active site of human AChE with a dimension of 40 × 30 × 40 Å (X, Y and Z axes of −13.987, −41.668 and 27.109, respectively). Molecular docking was consequently performed with AutoDock Vina (Trott & Olson, 2010) using default parameters. The docking protocol was validated in order to ensure its reliability for subsequent analysis of the studied compounds. This was performed by extracting the co-crystal ligand, donepezil, from the PDB file and re-docked to the co-crystal human AChE protein. The root mean squared deviation (RMSD) of the atomic position between the original orientation of the co-crystal ligand and the re-docked ligand is computed and is deemed acceptable if the RMSD value is less than or equal to 2.0 Å.

A set of 30 representative and chemically diverse compounds, which will be referred hereafter as the diversity set, were extracted from the full set of active AChE inhibitors (i.e., IC₅₀ <1 µM) using the Kennard–Stone algorithm (Kennard & Stone, 1969). These compounds were used as ligands for molecular docking against the human AChE. The binding energy (kcal/mol) of AChE inhibitors were calculated according to the built-in scoring function of Autodock Vina and conformers providing the lowest binding energy were selected for further analysis of the binding mode. Furthermore, key-interacting residues and their moiety preferences were analyzed using LigPlot+ (Wallace, Laskowski & Thornton, 1995), Maestro (Schrödinger, 2015b) and the SiMMap web server (Bollback, 2006). Finally, three-dimensional structure of protein–ligand interaction was created and visualized using Pymol (Schrödinger, 2015a).

Chemical space of AChE inhibitors

Navigation of the chemical space of AChE inhibitors was performed to gain insights into the structure–activity relationship by analyzing the Lipinski’s rule-of-five descriptors. Chemical space analysis may provide important knowledge on the general character of compounds governing inhibitory properties of compounds. Exploratory data analysis was performed using the Lipinski’s rule-of-five descriptors comprising of MW, ALogP, nHBDon and nHBAcc. MW represents the molecular size of a compound that is commonly used because of it can be easily interpreted and calculated as well as appropriate size of a compound is important for its passage via lipid membrane. ALogP is a widely used parameter for determining the lipophilicity of a compound and used in calculating the membrane penetration and permeability of compounds. nHBDon and nHBAcc describe the number of hydrogen bond donors and hydrogen bond acceptors, respectively, which is used to measuring hydrogen bonding capacity. Visualization of the chemical space of ALogP as a function of MW is shown in Fig. 2, as to investigate the chemical space of AChE inhibitors. A dense distribution of inhibitors was observed within the space of MW starting from approximately 300–600 Da and within the space of ALogP ranging from approximately −2.5 to 2.5. In addition, the box plot of the Lipinski’s descriptors is shown in Fig. 3. Compounds with negative ALogP values approximately of closer to 0.0 can be found in inactive inhibitors whereas most of the active inhibitors tend to possess approximately lower values in average of ALogP values.

Visual representation of the overall distribution of data values of Lipinski’s descriptors is shown as box plots in Fig. 3 in which the ALogP, nHBAcc, nHBDon and MW are shown in Figs. 3A, 3B, 3C and 3D, respectively. Analysis of the box plots revealed that there were no differences amongst the three bioactivity classes for nHBAcc and nHBAcc as deduced from the boundaries of the boxes (i.e., representing the first and third quartiles). ALogP and MW were found to display differences amongst the bioactivity classes. Particularly, ALogP values for actives were the lowest while negligible differences were observed for the other two classes. Furthermore, MW for actives were the largest amongst the three bioactivity classes, which is followed by the intermediates while inactives were smallest.

QSAR model for predicting AChE inhibitory activity

A data set comprising of 2,570 compounds were used for construction of QSAR models. Particularly, twelve sets of fingerprint descriptors were benchmarked in order to find the best performing set. Prior to modeling, feature selection was applied to remove collinear descriptors. Each of the twelve models were then built using a data split ratio of 80/20 in which 80% of the data set was used as the internal set and 20% as the external set. This procedure was iteratively performed in which each of the 100 independent data splits were used for model construction and the performance results given in Table 2 are the mean and standard deviation values derived from these runs.

It can be observed that all twelve models are capable of capturing the inhibitory activity space of AChE inhibitors as they provided R² and Q² (i.e., both 10-fold CV and external sets) greater than the threshold values proposed by Golbraikh & Tropsha (2002) of 0.6 and 0.5, respectively, which is indicative of robust model performance. The possibility of chance correlation can be assessed from the R²–Q² margin as described by Eriksson & Johansson (1996) in which values <0.2–0.3 are indicative of predictive and reliable models while values >0.2–0.3 suggests possible chance correlation or the presence of outliers in the data set. Furthermore, observation of the $Q_{\mathrm{CV}}^2$–$Q_{\mathrm{Ext}}^2$ margin revealed that the difference was negligible with values in the range of 0 and 0.01.

Generally, it can be seen that models with larger descriptor size, namely CDK and CDK extended, afforded the best performance with $Q_{\mathrm{CV}}^2$ of 0.79 ± 0.07 and 0.79 ± 0.06, respectively, and $Q_{\mathrm{Ext}}^2$ of 0.80 ± 0.04 and 0.81 ± 0.04, respectively. The opposite also holds true as the model with the least number of descriptors were also found to perform the worst amongst the other fingerprints with $Q_{\mathrm{CV}}^2$ of 0.55 ± 0.09 and $Q_{\mathrm{Ext}}^2$ of 0.56 ± 0.05. In a nutshell, the model performance in order of decreasing value is as follows: CDK extended > CDK > MACCS ≈ Substructure count ≈ Klekota–Rota count > PubChem > Klekota–Roth ≈ 2D atom pairs count > CDK graph only > Substructure > 2D atom pairs > E-state.

The best performing model is not necessarily the best choice considering the fact that the descriptor size for the best models were quite high and is consequently prone to overfitting. It was found that the substructure count provided reasonably good predictive performance (i.e., $Q_{\mathrm{CV}}^2$ and $Q_{\mathrm{Ext}}^2$ of 0.78 ± 0.06 and 0.78 ± 0.05, respectively) with the advantage of making use of a small set of 26 descriptors. Therefore, this fingerprint was selected for further interpretation of the feature importance.

To further check the reliability and validity of the selected model, Y-scrambling test was performed for 100 iterations. Table 3 demonstrates that QSAR models built using substructure count has a low Q² (−0.0013), which rules out the possibility of chance correlation. Furthermore, model afforded an $r_m^2$ value of 0.61 ± 0.06 thereby revealing its robustness. It is observed in Table 3 that the value of $\mathrm{\Delta }{r}_m^2$ is greater that 0.2 but also close to 0.20.

As shown in Fig. 4, it can also be seen that scatter plots of experimental versus predicted pIC₅₀ of panels A, C and F displayed narrower variance of the data points than the other methods as assessed via 10-fold cross-validation and external set.

Applicability domain

The AD of the proposed QSAR model was defined as provided by the Williams plot shown in Fig. 5. The employed data set consisting of 2,570 compounds was randomly split to two separate subset in which the first subset constituting 80% of the data set was used as an internal set while the second subset constituting the remaining 20% were used as an external set. Compounds representing the internal set (blue dots) and external set (red dots) are shown in the Williams plot and it can be clearly seen that almost all of the 2,570 compounds were located within the boundaries of applicability domain, which indicated that our proposed QSAR model had a well-defined AD.

As can be seen in Fig. 5, very few compounds indeed fall outside the ±3 standardized residual range. This consisted of six compounds (997, 1829, 62, 1096, 13, 677) from the internal set and 11 compounds (2116, 2120, 2388, 2117, 2392, 2323, 2423, 2424, 2507, 2219, 2422) from the external set that had standardized residual higher than 3. On the other hand, 7 compounds (1567, 576, 1644, 1098, 2022, 1447 and 322) from the internal set and 8 compounds (2486, 2353, 2130, 2553, 2389, 2072, 2103, 2125) from the external set had standardized residual lower than −3. The corresponding chemical structures are provided in Table S1.

Mechanistic interpretation of feature importance

Feature importance analysis help reveal features that are important toward bioactivity. There are essentially two parameters for evaluating the relative importance of features used in models using the RF algorithm: (i) accuracy and (ii) Gini index (i.e., variance of the responses). The latter was selected as a metric for ranking important features (i.e., mean decrease of the Gini index) for predicting the pIC₅₀ of AChE inhibitors (Fig. 6). Table 4 lists the substructure fingerprints along with their respective descriptions.

As can be seen in Fig. 6, the top ranking feature is secondary carbon (SubFPC2), which is a carbon atom with two carbon neighbors. In the context of drug design, such central carbon atom may be more difficult to be accessed and metabolized by cytochrome P450 (Uetrecht & Trager, 2007) and therefore are more metabolically stable.

The second most important feature is the rotatable bond (SubFPC302). Based on the rule of three for defining lead-like compounds, a compound may have a lead-like characters if it does have rotatable bonds of no more than 3. On the other hand, Veber et al. (2002) noted that the the upper limit of a orally bioavailable drugs is of seven rotatable bonds. Nevertheless, it has been found that number of rotatable bonds provide better discrimination between compounds that are orally active and those that are not. Kryger, Silman & Sussman (1999) claimed that E2020 (i.e., also known as donepezil and marketed as Aricept) needs at least two rotatable bonds on each side of the piperidine in which two aromatic moieties of E2020 interact with Trp86 and Trp286 (human AChE numbering), suggesting that links between aromatic systems of the inhibitor against its AChE counterparts are essential to yield high affinity.

The third important substructure is the aromatic ring (SubFPC274). Findings from X-ray crystallographic study showed that in the binding site of the co-crystal structure of AChE with tacrine, the aromatic ring of acridine engages in a π–π stacking interaction with the indole of Trp86 (human AChE numbering), thereby indicating the importance of the aromatic ring for AChE inhibition (Chen et al. 2012).

The fourth important feature is C ONS bond (SubFPC295), which is defined as the presence of any carbon connected with either oxygen, nitrogen or sulfur atom in a molecule. These atoms are considered as high electron density atoms, which exerted from higher electronegativity comparing with a carbon atom. Unequal sharing of electron pair making covalent bond contribute polarity and afford dipole moment to a molecule, which able to generate a dipole–dipole attraction such as hydrogen bond between two polar molecules. Furthermore, increasing of polarity and presenting of hydrogen bond improve water solubility, which is essential characteristic of a drug.

The fifth important substructure is secondary mixed amine (SubFPC32). The importance of the moiety was demonstrated in the work of Bembenek et al. (2008) in which a structure-based approach was used to reveal that in the catalytic triad, Trp86 interacts with the quaternary amine of ACh through a cation–π interaction. Furthermore, in the ‘anionic’ site Trp286 appears to attract the amine moiety via cation and/or hydrophobic interactions.

The sixth, seventh and eighth important substructures are heterocyclic rings (SubFPC275), tertiary carbon (SubFPC3) and primary carbon (SubFPC1). Heterocycles are of high relevance in the design of AChE inhibitors as it allows π–π stacking interaction with key amino acid residues in the binding site of AChE. It is observed that the binding site of the AChE are highly hydrophobic in nature. Particularly, the Trp286 (human AChE numbering) which is a part of the peripheral anionic site of the AChE is involved in the π–π interaction with heterocycles of AChE inhibitors (Lu et al., 2011). The aforementioned explanation made for the secondary carbon is also applicable for the tertiary carbon in which the higher number of carbon neighbors would also confer high stability against cytochrome P450 metabolism.

The ninth and tenth important substructures are Hetero N nonbasic (SubFPC181) and tertiary aliphatic amine (SubFPC26) in which the former is defined as aromatic nitrogen. The presence of the nitrogen atom, which is a cationic moiety deemed to interact with aromatic residues through π-cation interaction as observed in E2020 against Torpedo californica AChE (tcAChE). The charge of nitrogen atom located in piperidine ring provide π-cation interaction with the side chain of Phe337 (Phe330 in tcAChE numbering) (Guo et al., 2004). Since the active site gorge of AChE is comprised of several aromatic residues known as the aromatic patch, the addition of cationic moiety could possibly increase the binding affinity when the ligand is arranged in a suitable conformation with respect to the aromatic side chain of residues in the active site.

It is worthy to note that substructures pertaining to the covalent inhibitors, carbamates and organophosphates, were not found in the top ten important features. A manual inspection of the 307 substructure fingerprints revealed that there were none specifically describing the characteristic feature of carbamate and organophosphates. However, there were quite a few substructures that resembled partial features of the carbamate (e.g., carboxylic acid (SubFPC84), carboxylic ester (SubFPC85), “NOS methylen ester and similar” (SubFPC65), etc.) as well as substructures resembling partial features of organophosphates (e.g., phosphonic monoester (SubFPC230), phosphonic diester (SubFPC231), phosphonic monoamide (SubFPC232), phosphonic esteramide (SubFPC234), phosphonic acid derivative (SubFPC235), etc.). In spite of the presence of these descriptors, important features obtained from predictive model as revealed by the Gini index did not contain these features in the top ten. A possible explanation for such observation could be that there were a few carbamates (123 compounds) and organophosphates (18 compounds) present in the data set and its relative importance may have been masked by other features that represented the other structural class. Thus, it seems very lucrative for a future large-scale study to be performed focusing on these two compound classes.