2021 Volume 12 Issue 2

Theoretical Investigation of Some Donepezil-based Derivatives as Dual Inhibitors for beta-Amyloid- and Cholinesterase Enzymes

 

Assia Meziane, Amina Ghomri, Salim Bouchentouf*, Mohamed El-Shazly


Abstract

To treat Alzheimer’s Disease (AD), which is the most prevalent form of dementia, cholinesterase enzymes (AChE and BuChE) and amyloid-beta (Aβ) are attractive targets. In this work, different computational approach namely Density Functional Theory (DFT), Molecular Docking, and multi-QSAR modeling were performed on 22 donepezil-based derivatives which were reported as potent dual Aβ and (AChE and BuChE) inhibitors. The molecular geometries of the studied derivatives were carried out using GAUSSIAN 09 software with the level of theory (DFT, 6/31g*). The dual inhibitors adopted minimum energy. The results pointed out the importance of the inhibitors' geometries in enzyme inhibition. The QSAR models elaborated by means of Molecular Operating Environment (MOE) package, showed good statistical values for targets AChE (R²adj = 0.976, q2 = 0.871, RMS = 0.130), BuChE (R²adj = 0.976, q2 = 0.554, RMS = 0.092) and Aβ (R²adj = 0.861, q2 = 0.525, RMS = 0.113). To identify the binding pattern between the ligands and target enzymes, we implemented molecular docking studies for the datasets. The obtained information was related to the essential structural features that were related to the QSAR of the predicted models.

Keywords: Molecular docking, ButylCholinesterase, Quantitative structures activity relatioships, Acetylcholinesterase, Density functional theory, Donepezil


Introduction

The most common type of dementia is Alzheimer's Disease (AD) (Ronson, 2011; Uddin & Amran, 2018). AD affects over 44 million people worldwide. By 2030, the number of patients will be doubled and even tripled by 2050 (Harris, 2019). The worldwide cost of dementia was estimated at 604 billion USD in 2010 and is expected to increase exponentially. To tackle this problem, it is advisable to develop early detection methods and effective therapeutics (Tappen, 1997; Bamidis et al., 2020; Lilford & Hughes, 2020). AD may also affect patients with noncognitive disorders such as depression anxiety, hallucinations, and delusions (Vicente et al., 2015). Although AD pathogenesis is complex and unclear, there are several developed theories, but none of them revealed the specific cause of AD (Hardiman et al., 2016; Lilford & Hughes, 2020). Common targets were identified and explored over the last two decades. The main hallmarks such as amyloid plaques (Lee, 2000; Ronson, 2011), peptide aggregates, the tau protein aggregates (Stonebrook, 2007; Sigurdsson et al., 2012) discovered in AD patient brains were correlated with most of these targets. AD development was affected by the role of cholinergic deficit (Banner & Nixon, 1992; Piguet & Poindron, 2012). Aβ peptides are APP (Amyloid Precursor Protein) proteolytic by-products that consist of 42 (Aβ1−42) amino acids and40 (Aβ1−40). One of the strategies to cure AD is to block the generation of Aβ peptides aggregation. AChE, BuChE, and Aβ aggregation inhibitors emerged as effective tools for AD treatment. It was suggested that the effectiveness of the treatment was significantly improved by the dual inhibition strategy of these enzymes (Lajtha & Banik, 2001; Kuncharoenwirat et al., 2020; Sargazi & Taghian, 2020). The evaluations of in vitro Aβ aggregation inhibitors, progress, are time-consuming and labor extensive task. Computer-Aided Drug Design (CADD) (Kapetanovic, 2008) tools including molecular modeling in combination with Quantitative Structure-Activity Relationship (QSAR) (Li et al., 2019; Kasmi et al., 2020; Mahmud et al., 2020) and molecular docking (Roy et al., 2020) are used to evaluate a ligand activity and the kind of interactions into the protein active site (Taha & AlDamen, 2005). They provide useful tools for the design of new drugs to save time and money (Gu et al., 2021).

The biological activity of donepezil as anti-Alzheimer diseases has been the subject of different investigations. Shamsi et al. (2020), with human transferrin, applied molecular docking, calorimetric, and spectroscopic insights into the role of donepezil as an anti-Alzheimer’s drug. Aranda-Abreu et al. (2011) studied and discussed the absorption of donepezil after oral administration and compared it to tacrine. They found that donepezil was better tolerated by patients and caused fewer adverse reactions. Wallin et al. (2007) evaluated the results of three-year donepezil treatment and from their observation, the obtained results indicated a positive outcome in the routine clinical setting. To shed more light on donepezil derivatives as anti-Alzheimer’s drugs, we selected 22 derivatives (Table 1 and Figure 1) (Khosravan et al., 2017) and evaluated their cholinesterase enzymes (AChE and BuChE) inhibitory activity and amyloid-beta (Aβ) in silico (Yerdelen et al., 2015).

Materials and Methods

Dataset and Target Preparation

The molecular optimized geometries of the 22 donepezil derivatives were optimized using the DFT/ B3LYP (Daramola et al., 2010). Functional hybrid, with the 6-31g* basis set was implemented in the Gaussian 09 software (Hiscocks & Frisch, 2009), the stability of geometries was checked by the absence of the imaginary frequencies. The ligands properties (Table 1) were obtained by MOE software (Höltje, 2008; Kukol, 2008; Royal Society of Chemistry (Great Britain), 2014) results showed that the studied ligands were nontoxic, and the molecular weights were less than 500. The toxicity was done by MOE software.

The X-ray crystal structures of the targets, AChE (PDB ID: 1HBJ), BuChE (PDB ID: 4BDS), and Aβ (2BEG) were downloaded from the RCSB Database (http://www.rcsb.org/pdb). These three PDBs were chosen for modeling research since their crystal structures are in a condition that shows the pharmacological target for developing new drugs to treat AD.

QSAR Modeling and Models Validations

Biological Activities

Acetylcholinesterase (AChE), butyrylcholinesterase (BuChE), and Amyloid-beta (Aβ) inhibitory activities of a series of 22 donepezil-like amide secondary derivatives were taken from Kadir et al.’s work (Yerdelen et al., 2015), each activity was mentioned as IC50(μM) for (AChE) and (BuChE), and as a percentage for Aβ inhibitory activity. Values were expressed as mean ± standard error of the mean of three independent experiments. Values were converted to pIC50 as pIC50= −log IC50. The dataset was divided into a training set containing 17 compounds and a test set composed of 5 compounds (Table 1).

Molecular Descriptors Generation

Molecular descriptors were generated using MOE programs to predict the correlation between these parameters and their activities by developing a linear model (Partial Least Squares regression (PLS) (Wehrens & Mevik, 2007; Kovačević et al., 2018). In this work, 12 QSAR models were developed (Tables 2) using the experimental IC50 data. Calculations were done using a total of 365 different descriptors exploited in the MOE software. These sets of descriptors were first pre-processed with a variance threshold of 0.0001 and passed through a correlation coefficient of 0.99 to eliminate correlations between the input descriptors and noise level. A Genetic Algorithm (GA) was applied to select the best possible set of descriptors for QSAR modeling from the pre-processed pool of descriptors. This "descriptor elucidation" procedure allowed us to select the four most significant descriptors to build our QSAR models which were: AChE (chi0v_C, PEOE_VSA_POL, vsurf_D5, vsurf_Wp4), Bthe (npr1, vsurf_CP, vsurf_CW4, vsurf_Wp6) and for Aβ (a_ICM, density, vsurf_HL1, vsurf_ID1). The identified descriptors were the most pertinent for our elucidation as they reflect the required activities of all the studied molecules, assuming that a change of the molecular structure modifies the inhibitory activity of donepezil derivatives.

Regression Analysis

Data fitting was accomplished using PLS regression analysis. When there are many independent variables in the trial descriptor pool relative to the number of the dependent variables, the pIC50 endpoints, this data fitting technique becomes useful. When there is no method for ranking the individual members (molecular descriptors) of the trial descriptor pool and/or knowing possible inter-relationships among the training set of descriptors, the PLS is applied. Regression analysis was done by PLS method, because of the largest of the independent variables in the trial descriptor compared with the number of dependent variables, pIC50, our choice is justified because PLS is the best one as there is no method for ranking the individual members (descriptors) or distinguishing the inter-relationships between those descriptors. The PLS regression was used to build a single QSAR model that contains the entire trial descriptor pool. Both R2 and the leave-one-out (LOO) cross-validated correlation coefficient, q2, were employed to characterize the quality of the resultant QSAR models.

Molecular Docking Procedure

Using MOE, molecular docking studies were done aiming to find the best conformation of the donepezil derivatives in protein binding sites. The downloaded 3D structures of the protein complexes were protonated, and energy minimized in an MMFF94x force field to a gradient of 0.0001 kcal/mol/Å (Engh & Huber, 1991). Using the MOE-Alpha Site finder, the active sites were generated. The dummy atoms were created from the obtained alpha spheres. Prepared protein structures are presented in Figure 2 for AChE, Bthe and Aβ prepared enzymes and the active sites are also presented. The 3D structures of all compounds were built by GaussView06, structures were optimized and saved as Mol Folders. Using the MOE program, the database set was created, and this database was used as input MOE-docking. We used the default settings of the parameters with Ligand Placement (Triangle Matcher) and Rescoring (London dG) implemented in the MOE program (Chemical Computing Group Inc, 2016). Furthermore, the LigX feature of MOE was used to find the hydrogen bonding interactions between the ligand and receptor protein. The best docking scores were used for the calculation of binding energy.

Results and Discussion

DFT Calculations and QSAR Modeling

To build the QSAR models, molecular geometries were optimized using the B3LYP/6-31g* and the database was created. The elaborated QSAR models (Tables) revealing the correlations between the AChE, BChE, and Aβ inhibitory activities, and the corresponding equations are also presented in the above tables.

The terms were defined including NTraining the number of molecules in the training set, NTest the number of molecules in the test set, R²CV (Q2) the LOO cross-validated coefficient, RMSE the root mean square error, and R² the correlation coefficient. The absolute difference between the activity field and the value of the model is $Z-SCORE, which is divided by the square root of the mean square error of the data set. For the external validation, the values of R²test (correlation coefficient), RMSE test (Root Mean Square Error) were selected corresponding to the best models for each enzyme.

R2 allowed us to compare the experimental and predicted studied activities, which is the protein inhibition. A good model must have a value of R2 above or equal to 0.5 whereas RMSE is mostly used to decide if the QSAR model possesses the predictive quality reflected in R2. It showed the error between the mean of the experimental values and the predicted properties. If RMSE is above 1 (RMSE ˃1), the model has a poor ability to foretell the properties even with a good R2 value. However, R2 and RMSE parameters are not sufficient to judge the QSAR validity, that is why cross-validation is required. Cross-validation R2cv is used to allow the determination of how large the model can be used for a random data set and to evaluate the predictive power of the model’s equation. On the other hand, the standardized value, which specifies the exact location of an X value within a distribution by describing its distance from the mean in terms of the standard deviation units is represented by XZSCORE and ZSCORE. This subset must be examined carefully to detect errors or to determine new descriptors to be calculated. ZSCORE must be less than 2.5 to suggest that the QSAR model is good to be used. Internal validation is not sufficient to estimate the predictive power of a QSAR mode.l. Golbraikh and Tropsha suggested (Golbraikh et al., 2003) the following statistical characteristics of the test set: correlation coefficient R between the predicted and observed activities, coefficients of determination (R2test) the correlation between the predicted and the experimental inhibition activities are represented in Figure 3 for AChE, BuChE, and Aβ respectively. They consider a QSAR valid model with R2test > 0.6 taken as an indicator of good external predictability.

To highlight the validity of our models for each enzyme as well as the chosen descriptors contributions on each activity, we calculated the statistical parameters shown in (Tables 2 and Figures 3) for the internal validation and correlation coefficients for the external validation (Figure 3).

In the case of AChE, Table 1 and Figure 3 showed that among the four elaborated models, model number 2 demonstrated the best results with the value R2 = 0.88 and RMSE = 0.13. For the cross-validation and R2CV = 0.58 with $Z-SCORE less than 2.5 (0.6946 < ZS < 1.142), our results allowed us to confirm the correlation between the six selected descriptors (E_str; PM3_IP; PM3_LUMO; SMR_VSA2; SMR_VSA3; vsurf_Wp4) and the donepezil derivatives. AChE inhibition specified as IC50 correlation plot between the predicted and experimental activities for AChE enzyme confirmed that model 2 is the best. Since the internal validation is insufficient to judge the model validity, we proceeded to an external validation using five molecules from the dataset which were not used in the model elaboration. Figure 3 gave the values and correlations between the experimental and theoretical activities obtained for the five test molecules. Our results showed that the selected model 2 can be used to predict the donepezil derivatives activities. The residual values between our predicted values and the experimental ones varied from 0.0013 to 0.1189. Those values are considered excellent values for the prediction of inhibitory activities. As indicated by Jalali-Heravi and Kyani, we noted that the positive and negative values varying the residuals on both sides of zero showed that there was no systemic error (Jalali-Heravi & Kyani, 2004). To predict the pAChE inhibitory activities, this model can be applied successfully. The correlation coefficient between the predicted and experimental values is R2 = 0.96 and this value was also considered excellent.

The selected descriptors chosen in the best model 2 suggested that E_str, PM3_LUMO, and vsurf_Wp4 were associated with a negative coefficient, indicating that increasing the value of these descriptors was unfavorable to the acetylcholinesterase inhibitory activity for each molecule PM3_IP, SMR_VSA2, and SMR_VSA3 descriptors that were associated with a positive coefficient suggesting that the inhibitory activity of the molecules increased with the increase of the approximate accessible van der Waals surface area of the molecule. The same conclusions could be made for the ionization potential of donepezil derivatives.

In case of BuChE, Table 2 and plot 2 showed that among the four elaborated models, model number 6 is best model showing the best results with values of R2 = 0.74, RMSE = 0.09 and the cross-validation correlation coefficient R2CV = 0.55 with $Z-SCORE less than 2.5 (0.52 < ZS < 0.71). Between the six selected descriptors. E_oop, MNDO_dipole, vsurf_DD12, and the vsurf_ID4 were associated with negative coefficients, so an increase in these quantities decreased the inhibitory activities. The increase of the vsurf_HB3 and vsurf_ID3 corresponded to an increase of the inhibitory activity. These descriptors are useful in pharmacokinetic property prediction. The correlation plot between the predicted and experimental activities for the BuChE enzyme is presented in Figure 3 and it confirms that model 6 is the best one between the four established models. For the external validation, we used five molecules from the dataset which were not used in the model elaboration. Table 2 and plot 2 provided the values and correlations between the experimental and theoretical activities obtained for the five tested molecules. Our outcomes revealed that the selected model 6 can be used to predict the donepezil derivatives inhibitory activities. The residual value between our predicted values and the experimental ones varied from 0.039 to 0.06. These values are also considered excellent values for the prediction of BuChE inhibitory activity. The correlation coefficient between the predicted and experimental values (plot) was R2test = 0.95, which was also classified as an excellent value.

In the case of Aβ, Table 2 and Figure 3 showed that among the four elaborated models, model number 11 is the best model with R2 = 0.76 and RMSE = 0.11. For the cross-validation and R2CV = 0.52 with $Z-SCORE less than 2.5 (0.01 < ZS < 2.3), our results allowed us to confirm the good correlation between the six selected descriptors (a_ICM; density; vsurf_HL1; vsurf_ID1; vsurf_ID2) and the Aβ inhibition by donepezil derivatives specified as IC50. All parameters were expressed with positive coefficients, suggesting that an increase in the values of these descriptors will lead to an increase in the inhibitory activity. The correlation plot between the predicted and experimental activities for the Aβ enzyme is presented in Figure 3 and it confirms that model 11 is the best. We also proceeded to an external validation using five molecules from the dataset which were not used in the model development. Table 2 and Figure 3 provided the values and correlations between the experimental and theoretical activities obtained for the five test molecules. Our results suggested that the selected model 11 can be used to predict the donepezil derivatives activities. The predicted values were very close to the experimental values and they were considered excellent values. The correlation coefficient between the predicted and experimental values (Figure 3) was R2 = 0.90, which was also considered an excellent value.

Molecular Docking Simulation

To prepare the enzymes, we identified the active site residues of both cholinergic and amyloid-beta targets, the active-site residues in AChE (Figures 4) are GLY118, GLY119, and ALA 201, which create the catalytic triad. The active-site residues for BuChE (Figures 4): HIS438, SER198, and GLU325. For the Aβ (Figures 4) the active site did not show the potent residues because the downloaded structure was not complexed with any reference ligand. All the prepared enzyme structures and the fixed sites are presented in Figures 4.

Studies on molecular docking proceeded with the 22 molecules database designed by MOE. The AChE, BuChE, and Aβ protein-ligand complexes were evaluated by estimating various types of interactions (non-polar and polar such as H-bonding interactions, electrostatic interactions, van der Waal’s interactions, hydrophobic interactions for ligands). The results indicated that the lowest docking score AChE (-8.1299 kcal/mol, ligand 4), BuChE (-7.0399 kcal/mol, ligand 14), and Aβ (-4.8651 kcal/mol, ligand 15) were the best docking score and other scores were calcified respecting this order (that all scores were calcified and the lowest founded for each enzyme was the one given). Figure 4 showed the 2D Protein-ligand interaction maps for the 22 molecules. The figure showed that the best dual ligands that gave good scores with the three studied targets and can be used as a model to design new dual inhibitors are ligand number 19 for AChE and BuChE duality and ligand number 15 for the BuChE and Aβ duality.

 

 

Figure 1. The Molecular Structures Corresponding to the 22 Ligands

Table 1. Molecules Names and IC50 against AChE, BthE, and Aβ and obtained Ligands Properties using MOE of Compounds (1-22)

 

Model Equation

 

Ache pic50 =

1.17613

-0.14325 * chi0v_C

+0.06162 * PEOE_VSA_POL

-0.16995 * vsurf_Wp4

Ache pic50 =

-1.18719

-0.15561 * E_str

+0.29230 * PM3_IP

-0.17014 * PM3_LUMO

+0.09599 * SMR_VSA2

+0.04150 * SMR_VSA3

-0.13785 * vsurf_Wp4

Ache pic50 =

1.56209

-0.13651 * E_str

-0.17214 * PM3_LUMO

+0.10080 * SMR_VSA2

+0.02918 * SMR_VSA3

-0.13451 * vsurf_Wp4

Ache pic50 =

-2.80053

+0.17094 * MNDO_dipole

+0.01555 * PEOE_VSA+0

+0.06022 * PEOE_VSA_HYD

+0.12946 * PEOE_VSA_POL

-0.02129 * PEOE_VSA_POS

-0.03965 * vdw_vol

 

 

BuChE

 

Validation

Metrics

Model 5

Model 6

Model 7

Model 8

Threshold

 

Internal

N Training

17

17

17

17

 

 

 

N Test

5

5

5

5

 

 

 

R2

0,69735

0,74745

0,70206

.71239

> 0.5

 

 

RMSE

0,09627

0,09220

0,09981

0,09170

 

 

 

R2cv

0.51263

0.55450

0.51057

0.53343

> 0.5

 

 

$Z-SCORE

0.4548<ZS<0.6322

0.7161<ZS<0.5254

0.4136<ZS<0.1081

0.5859<ZS<0.9151

< 2.5

 

External

R2test0.9566

 

RMSE test

 

0.9781

 

 

 

 

Model Equation

 

Bthe pic50 =

-1.95177

-0.02273 * MNDO_dipole

-0.04080 * vsurf_DD12

-0.01469 * vsurf_HB3

+0.15864 * vsurf_ID3

+0.15076 * vsurf_ID4

+0.01556 * vsurf_W3

Bthe pic50 =

-0.86927

-0.03012 * E_oop

-0.02579 * MNDO_dipole

-0.04784 * vsurf_DD12

+0.00154 * vsurf_HB3

+0.42363 * vsurf_ID3

-0.25670 * vsurf_ID4

Bthe pic50 =

-1.15305

-0.01916 * E_ele

+0.00508 * MNDO_dipole

-0.04552 * vsurf_DD12

+0.00060 * vsurf_HB3

+0.06793 * vsurf_ID3

+0.07737 * vsurf_ID4

Bthe pic50 =

-0.97705

-1.79169 * npr1

+1.09797 * vsurf_CP

+1.14182 * vsurf_CW4

-0.70721 * vsurf_Wp6

 

 

 

Validation

Metrics

Model 9

Model 10

Model 11

Model 12

Threshold

 

Internal

N Training

17

17

17

17

 

N Test

5

5

5

5

 

R2

0,75408

0,75952

0,76745

0,74128

>0.5

RMSE

0,11735

0,11516

0,11840

0,11386

 

R2cv

0.49967

0.51508

0.52146

0.52535

 

$Z-SCORE

2.10945<ZS<0.3466

2.0979<ZS<0.5309

1.3243<ZS<0.2875

2.3291<ZS<0.4602

<2.5

External

R2test 0.9011

RMSE test

 

 

 

0.9492

 

Model Equation

 

Aβpic50 =

-1.42650

-1.98824 * a_ICM

+1.88270 * density

+11.32901 * vsurf_HL1

+0.25480 * vsurf_ID1

+0.11959 * vsurf_ID6

Aβpic50 =

-1.42270

-2.13016 * a_ICM

+2.18384 * density

+11.71892 * vsurf_HL1

+0.68498 * vsurf_ID1

-0.30229 * vsurf_ID2

Aβpic50 =

-2.00729

-2.13629 * a_ICM

+2.50472 * density

+11.00672 * vsurf_HL1

+0.41933 * vsurf_ID1

+0.24070 * vsurf_ID8

Aβpic50 =

-1.49915

-2.08249 * a_ICM

+2.18037 * density

+11.48446 * vsurf_HL1

+0.40971 * vsurf_ID1

 

                         

 

 

a) Model 2 (AChE)

b) Model 6 (BuChE)

c) Model 11(Aβ)

Figure 3. Correlation Plots between Predicted and Experimental Activities Expressed pIC50for AChE, BuChE and Aβ

 

Ligands Interactions Maps between the Best Ligands Scored Ligands and BuChE

1Lig(14)

2Lig(22)

3Lig(15)

4 Lig(19)

5Lig(9)

6Lig(10)

Ligands Interactions Maps between the Best Ligands Scored Ligands and Aβ

1 lig(15)

2 lig(8)

3 lig(7)