Docking and Quantitative Structure Activity Relationship studies of Acyl Guanidines as β-Secretase (BACE1) Inhibitor

Abstract

${\beta}$-Secretase (beta-amyloid converting enzyme 1 [BACE1]) is involved in the first and rate-limiting step of ${\beta}$-amyloid ($A{\beta}$) peptides production, which leads to the pathogenesis of Alzheimer's disease(AD). Therefore, inhibition of BACE1 activity has become an efficient approach for the treatment of AD. Ligand-based and docking-based 3D-quantitative structure-activity relationship (3D-QSAR) studies of acyl guanidine analogues were performed with comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA) to obtain insights for designing novel potent BACE1 inhibitors. We obtained highly reliable and predictive CoMSIA models with a cross-validated $q^2$ value of 0.725 and a predictive coefficient $r{^2}_{pred}$ value of 0.956. CoMSIA contour maps showed the structural requirements for potent activity. 3D-QSAR analysis suggested that an acyl guanidine and an amide group in the $R_6$ substituent would be important moieties for potent activity. Moreover, the introduction of small hydrophobic groups in the phenyl ring and hydrogen bond donor groups in 3,5-dichlorophenyl ring could increase biological activity.

Introduction

Alzheimer’s disease (AD) is the most common neurodegenerative disorder and a major cause of dementia in the elderly; AD currently affects more than 30 million people worldwide.1,2 The number of AD patients is expected to increase significantly with aging of the human population. AD causes progressive damage to the nerve systems responsible for memory and is associated with the formation of insoluble amyloid plaques composed of β-amyloid (Aβ) peptides.3-5 Currently, the approved treatments for AD patients are acetylcholinesterase inhibitors and N-methyl-d-aspartate (NMDA) receptor antagonists. However, they only treat disease symptoms; disease progression is not treated.

The Aβ-peptides are generated from the sequential proteolysis of amyloid precursor protein (APP) first by β-secretase (beta-amyloid converting enzyme 1 [BACE1]), followed by γ-secretase.6-9 Because BACE1 is responsible for the first and rate-limiting step in the production of Aβ, it is an attractive therapeutic target for AD. It has been also reported that production of all forms of Aβ decreased in BACE1-knockout mice, although the mice were healthy.10,11 Therefore, BACE1 inhibitors have been considered as the preferred therapeutic agents for lowering brain Aβ levels in the treatment or prevention of AD.

There have been numerous reports on BACE1 inhibitors by many pharmaceutical companies, and some compounds have been tested in clinical trials.12-15 BACE1 is a membrane protein found primarily in the brain and is a member of the pepsin family under the aspartyl protease superfamily. BACE1 shares sequence homology with BACE2 (52%), cathepsin D (29%), pepsin (27%), cathepsin E (27%), and renin (24%).16 Therefore, BACE1 inhibitors should not only demonstrate potent activity but also selectivity, good oral bioavailability, blood–brain barrier penetration, and efficacy at a therapeutically acceptable dose. To overcome this challenge, 3D-QSAR approaches have been a strong research tool in medicinal chemistry and critical to induce some rational guidelines for further structural modification.

In this paper, we report ligand-based and docking-based 3D-QSAR analysis of acyl guanidine derivatives, which have been reported as potent BACE1 inhibitors.

 

Experimental

Data Set. The structures and biological activities data of the compounds used in this study were obtained from literature.17 Biological activity data were given as Ki values, which were converted to negative logarithmic scale (pKi), and were subsequently used as the dependent variables for the QSAR analyses. Fifty-two compounds were divided into a training set (40 compounds) and a test set (12 compounds). The 3D-QSAR models obtained from the training set were validated by the selected test set compounds. Structures and their biological activities are given in Tables 1, 2, and 3.

Table 1.aTest set compounds. This table is cited from reference 17.

Table 2.aTest set compounds. This table is cited from reference 17.

Table 3.aTest set compounds. This table is cited from reference 17.

Molecular Structures. All computational studies were performed with the SYBYL-X 1.3 molecular modeling software package.18 The structures of compounds were generated with a sketch tool and energy minimization was carried out using a TRIPOS force field with the Powell method and conjugate gradient termination. Partial atomic charges of molecules were calculated based on the Gasteiger-Hückel charges. Low-energy conformation was searched by the simulated annealing method, and molecular alignment was achieved by the Distill rigid method in SYBYL. The most active compound (52) in the training set was used as a template molecule and an acyl guanidine moiety was used as a common substructure in the alignment.

Docking Studies. The crystal structure of BACE1 was obtained from the Protein Data Bank (PDB code: 4FSL), and docking analysis was performed using the Surflex-dock module to acquire the probable bioactive conformations of ligands.19 The PDB structure was downloaded and prepared using Preparation Tool in SYBYL, which was performed with the following steps: removal of waters, ligands, and cofactors; addition of hydrogen; and charging of protein with AMBER7 FF99. These steps were then followed by minimization of the protein with an AMBER7 FF99 force field. The compounds were docked into the protein using ligand mode with an empirical scoring function. Top 50 docked poses ranked by total scores were obtained for each compound. After top 5 docked poses were chosen from the 50 docked poses, each pose was put into data set one by one and PLS analysis was performed. The bioactive conformation for each compounds was selected from the pose which exhibit best results of QSAR models.

CoMFA and CoMSIA Models. CoMFA was carried out on the steric and electrostatic fields with the default values. 3D cubic lattices with a grid spacing of 2.0 Å were created automatically around the aligned molecules. The steric and electrostatic field energies of CoMFA were calculated for each molecule by using Lennard-Jones potential and Coulombic potential, respectively. A default sp3 carbon atom with a charge of +1 and a Van der Waals radius of 1.52 Å was used as a probe to generate the CoMFA steric and electrostatic fields. The computed field energy was truncated to 30 kcal/mol for both fields.

The CoMSIA model evaluates the following five physicochemical properties: steric, electrostatic, hydrophobic, hydrogen bond acceptor, and hydrogen bond donor fields. The CoMSIA method involves a common probe atom and similarity indices determined in regularly spaced grid points for the aligned molecules. The common probe atom with a radius of 1.0 Å, charge of +1, hydrophobicity of +1, hydrogen bond donating of +1, and hydrogen bond accepting of +1 was used to calculate the five fields. A default value of 0.3 was used for the attenuation factor.

Partial Least Squares (PLS) Analysis. A partial least squares (PLS) method was conducted with cross-validation to determine the optimum number of components, which were then used for the final 3D-QSAR models by non-crossvalidation.20 Cross-validation was performed with the leaveone-out (LOO) method in which one compound was removed from the data set and its biological activity was predicted with the model derived from the rest of the data set. Finally, non-cross-validated analysis was carried out using the optimal number of components for CoMFA and CoMSIA. To test the predictive power of 3D-QSAR model as a predictive tool, a test set with 12 compounds was used to predict biological activities.

PLS regression analysis was used to determine the linear correlation between the CoMFA and CoMSIA descriptors as independent variables and biological activities values as dependent variables. The predictivity of the models is evaluated by the r2pred value, which is expressed by the following function:

where SD is the sum of the squared deviations between the biological activities of the test set and the mean activity of the training set, and PRESS is the sum of squared deviations between the actual and predicted activity values of the test set.

 

Results and Discussion

3D-QSAR Analysis. CoMFA and CoMSIA analysis were derived for a data set of 52 acyl guanidine analogues as β-secretase (BACE1) inhibitors that were divided into a training set containing 40 compounds for model construction and a test set with the remaining 12 compounds for external model validation. The pKi values spanned about 3.4 orders of magnitude and were distributed across a range of values (Tables 1, 2, and 3). The compounds of the training and test sets were selected in order to cover the entire range of pKi values and structures.

Ligand-based CoMFA and CoMSIA Models. The statistical parameters of the CoMFA and CoMSIA models are summarized in Table 4. The correlation plots between the predicted and the actual pKi values of the training, and test sets are depicted in Figure 1(a) and 1(b). In the CoMFA model, PLS analysis showed a cross-validated q2 of 0.570 with four principal components and a non-cross-validated r2 value of 0.929. The F test value and standard error of estimation (SEE) were 114.509 and 0.285, respectively. The CoMSIA model with electrostatic, hydrophobic, and hydrogen bond donor fields gave a cross-validated q2 value of 0.690, a non-cross-validated r2 value of 0.936, and PLS components of 4. The F value was 128.308 and standard error of estimation (SEE) was 0.270.

Table 4.ONCa = optimum number of components; q2b = leave-one-out (LOO) cross-validated correlation coefficient; r2c = non-cross-validated correlation coefficient; SEEd = standard error of estimate; Fe = F test value; Sf = steric; Eg = electrostatic; Hh = hydrophobic; Di = donor; Aj = acceptor.

Figure 1.Correlation plots between the actual and predicted activities of 3D-QSAR models based on ligand alignment. (a) CoMFA model and (b) CoMSIA model. pKi values were converted into pKi (-log pKi) values (● : training set compounds, ▲: test set compounds).

Docking-based CoMFA and CoMSIA Models. The statistical results for the CoMFA and CoMSIA models are listed in Table 5. The CoMFA model with both steric and electrostatic fields gave a cross-validated correlation coefficient q2 value of 0.52 with an optimized component of 4. This model obtained a non-cross-validated correlation coefficient r2 value of 0.952, an F value of 172.177, and a standard error estimate of 0.239. The contributions of steric and electrostatic fields were 44.0% and 56.0%, respectively.

Table 5.ONCa = optimum number of components; q2b = leave-one-out (LOO) cross-validated correlation coefficient; r2c = non-cross-validated correlation coefficient; SEEd = standard error of estimate; Fe = F test value; r2predf= predictive correlation coefficient; Sg = steric; Eh = electrostatic; Hi= hydrophobic; Dj = donor; Ak = acceptor.

The CoMSIA model using electrostatic, hydrophobic, and hydrogen bond donor fields provided a good cross-validated correlation coefficient q2 of 0.725 with an optimal component number of 4, which indicated that the model could be a useful tool for predicting the pKi values. A high non-crossvalidated correlation coefficient r2 of 0.958 with a low standard error estimate of 0.222 and an F value of 200.981 were obtained. The corresponding contributions of electrostatic, hydrophobic, hydrogen bond donor fields were 38.1%, 29.6%, and 32.3%, respectively. Therefore, the electrostatic field had a strong influence on the biological activity, and the hydrophobic and hydrogen bond donor fields showed almost the same influence to the activity. These results suggested that the electrostatic interactions of the molecules with the receptor could be the most important for the BACE1 inhibiting activity; the hydrophobic and hydrogen-donating interactions were also important.

As shown in Tables 4 and 5, the CoMSIA model based on docking alignment exhibited better statistical parameters than those from other models. External validation was carried out to evaluate the external predictivity of the docking-based 3D-QSAR models. Twelve acyl guanidine analogues, which were not included in the training set, were used to assess the predictive ability of the CoMFA and CoMSIA models. The predictive correlation coefficients r2pred were 0.929 and 0.956 for the CoMFA and CoMSIA models, respectively, which suggested that the CoMSIA model based on the docking method was reliable and able to predict binding affinities of new analogues accurately.

From this CoMSIA model, the predicted pKi values and the residuals of the training set and test set were calculated and presented in Table 6. The graphs of actual versus predicted pKi values for the training and test sets are shown in Figure 2(a) and 2(b). The predicted pKi values were well correlated with actual values and showed small residual values, which suggested the good predictive ability of the model.

Table 6.The actual pKi, predicted pKi, and residuals of the training set and test set compounds based on docking alignment in the CoMSIA model

Figure 2.Correlation plots between the actual and predicted activities of 3D-QSAR models based on docking alignment. (a) CoMFA model and (b) CoMSIA model. pKi values were transformed into pKi (-log pKi) values (●: training set compounds, ▲: test set compounds).

3D-contour maps were used to visualize information of the derived 3D-QSAR models (Fig. 3). The maps use the characteristics of the compounds that are important for biological activity and show the regions around molecules where increased or decreased activities are expected based on physicochemical property changes in the molecules. In an electrostatic contour map, the blue and red regions illustrate the favorable sites for electropositive and electronegative groups, respectively. The yellow contours indicate where hydrophobic groups enhance the biological activity, whereas the white contours show regions where hydrophilic groups increase activity.

Figure 3.3D-contour maps of the CoMSIA models with electrostatic, hydrophobic, and hydrogen bond donor fields. (a) Contour map of electrostatic field; (b) contour map of hydrophobic field; (c) contour map of hydrogen-bond donor field. Compound (51) is shown within the fields (blue, favorable positive charge; red, favorable negative charge; yellow, favorable hydrophobicity; white, unfavorable hydrophobicity; purple, favorable hydrogen bond donor; cyan, unfavorable hydrogen bond donor).

The purple area favors hydrogen bond donor groups, while the cyan represents an unfavorable region for hydrogen bond donor groups. The molecule in the contour maps is compound (51) which showed the most potent activity.

The CoMSIA electrostatic contour map is shown in Figure 3(a). Blue contours were observed near the NH of the amide group in the R6 substituent, the 5-methyl group in the isothiazole ring, and below the NH of the guanidine. This indicated that positively charged groups in these areas would increase the activities. The red contours near the carbonyl of the acyl guanidine group and the CH2 next to the amide group in R6 showed that negatively charged groups in these regions were favorable for the activity. As shown in Figure 3(a), compound (51) had blue contours around the NH and red contours near the carbonyl group. These contours were in accordance with the electrostatic contour maps to show potent activity for compound (51). In compounds (32 and 47), red contours were in oxygen atoms of carbonyl group and blue contours were found around nitrogen atoms of amid group to give potent activities. However, compounds (8, 11, and 28) were not corresponded to the electrostatic contour map to present the low biological activities, in which they had blue contours near negatively charged groups (F, Cl, and carbonyl group).

The CoMSIA hydrophobic contours are illustrated in Figure 3(b). There were yellow contours around the benzene ring and two chlorides where hydrophobic interactions enhanced the biological activity. The white contours were close to the amide group of the R6 substituent and methoxy group where hydrophilic interactions increased the activity. In compounds (32, 51, and 52), yellow contours were near the hydrophobic groups (benzene ring and methyl group) and white contours were shown around the hydrophilic amide group of R6 substituent. These compounds were well correlated with the hydrophobic contour map to give increased activities. But the hydrophilic carbonyl groups were in yellow contour and hydrophobic benzene rings were close to white contour in compounds (4 and 5) to show very low activities.

The hydrogen bond donor contour map is shown in Figure 3(c). The cyan contour was present near the carbonyl of the amide group in R6 where hydrogen bond donor groups were unfavorable for the activity. The three purple contours were positioned around the NH and NH2 of the guanidine group and C2 position of 3,5-dichlorophenyl ring where hydrogen bond donor groups were favorable for activity. Compound (51) was in good correspondence with the cyan and purple contours to exhibit potent activity. The acyl guanidine group and an amide group in the R6 substituent were close to purple contours to give hydrogen bond binding site to carboxylic groups of the aspartates in the receptor.

The binding mode between compound (51) and receptor is depicted in Figure 4. The purple contours in Figure 3(c), were close to the carboxylic groups of the aspartates in the receptor (Fig. 4) to form hydrogen bonds between the NH and NH2 groups of guanidine and the carboxylic groups of the Asp80 and Asp276 residues.

Figure 4.The binding mode between compound 51 and BACE1

Another hydrogen bond was found between the NH of the amide group in R6 and the carbonyl group in the Phe156 residue. Compound (33) in X-ray crystal structure (4FLS) also showed a binding mode similar to that shown by compound (51), and this could explain why these compounds (51 and 33) are potent inhibitors of BACE1 (4 and 5 nM).

From the 3D-QSAR analysis, following structure modification was suggested to design the more potent BACE1 inhibitors (Fig. 5). It was suggested that the introduction of hydrophobic groups to C2 position of methoxy phenyl ring, hydrogen bond donor groups to C2 position of 3,5-dichlorophenyl ring, and negatively charged groups to methylene hydrogen of CH2 azetidinyl group might increase biological activity. It was also recommended that the substitution of positively charged groups for 5-methyl group in isothiazole ring and hydrophilic groups for methoxy group could be helpful for potent activity. However, an acyl guanidine group and an amide group in the R6 substituent were involved in hydrogen bonds between ligand and receptor. Therefore, these groups would be important moieties for the biological activity and their structural modification would not be recommended.

Figure 5.Design of new BACE1 inhibitors.

 

Conclusion

In this research, molecular docking studies were performed to explore the interaction mode between BACE1 and acyl guanidine inhibitors. 3D-QSAR studies helped to construct highly accurate and predictive 3D-QSAR models, including the CoMFA (q2, 0.520; r2, 0.952; r2 pred, 0.929) and CoMSIA (q2, 0.725; r2, 0.958; r2 pred, 0.956) models, based on docking alignment to predict the biological activity of new compounds. 3D-QSAR analysis suggested that acyl guanidine and amide groups in the R6 substituent were correlated well with contour maps and they would be important moieties for biological activity. It was suggested that the introduction of hydrophobic groups to C2 position of methoxy phenyl ring, hydrogen bond donor groups to C2 position of 3,5-dichlorophenyl ring and the substitution of positively charged groups for 5-methyl group in isothiazole ring, hydrophilic groups for methoxy group could be helpful for potent activity. Furthermore, the combined docking and 3D-QSAR studies may help to interpret the structure activity relationship of BACE1 inhibitors and to provide valuable insights into rational drug design for further improvement of the biological activity of BACE1 inhibitors.