Inhibition of acetylcholinesterase by two genistein derivatives: kinetic analysis, molecular docking and molecular dynamics simulation.

In this study two genistein derivatives (G1 and G2) are reported as inhibitors of acetylcholinesterase (AChE) and butyrylcholinesterase (BuChE), and differences in the inhibition of AChE are described. Although they differ in structure by a single methyl group, the inhibitory effect of G1 (IC50=264 nmol/L) on AChE was 80 times stronger than that of G2 (IC50=21,210 nmol/L). Enzyme-kinetic analysis, molecular docking and molecular dynamics (MD) simulations were conducted to better understand the molecular basis for this difference. The results obtained by kinetic analysis demonstrated that G1 can interact with both the catalytic active site and peripheral anionic site of AChE. The predicted binding free energies of two complexes calculated by the molecular mechanics/generalized born surface area (MM/GBSA) method were consistent with the experimental data. The analysis of the individual energy terms suggested that a difference between the net electrostatic contributions (ΔE ele+ΔG GB) was responsible for the binding affinities of these two inhibitors. Additionally, analysis of the molecular mechanics and MM/GBSA free energy decomposition revealed that the difference between G1 and G2 originated from interactions with Tyr124, Glu292, Val294 and Phe338 of AChE. In conclusion, the results reveal significant differences at the molecular level in the mechanism of inhibition of AChE by these structurally related compounds.

Introduction

Alzheimer׳s disease (AD), a progressively degenerative disorder of the brain, is believed to be a multifactorial syndrome with several target proteins contributing to its etiology. AD is characterized by a loss of basal forebrain neurons and reduced cortical and hippocampal levels of acetylcholine (ACh). The relation between the observed cholinergic dysfunction and AD severity provides a rationale for the therapeutic use of acetylcholinesterase inhibitors (AChEIs).

The inhibition of acetylcholinesterase (AChE, E.C. 3.1.1.7), which is responsible for the breakdown of ACh, has been proven as a successful way to relieve some cognitive and behavioral symptoms of AD3, 4. So far, several AChE inhibitors (Fig. 1), such as tacrine, galanthamine, huperzine A and donepezil have been used mainly for the clinical treatment of AD, all of which slow down neurodegeneration in AD patients to some extent.

Figure 1

The structures of AChE inhibitors.

The three-dimensional structure of AChE, as determined by X-ray crystallography for a large number of enzyme–ligand complexes9, 10, 11, reveals two main binding sites: the catalytic active site (CAS), comprising the Ser-His-Glu catalytic triad, and the peripheral anionic site (PAS), connected by a deep, hydrophobic gorge.

Tacrine is the first AChE inhibitor permitted by the FDA. The co-crystal structure of AChE (protein data bank entry:1ACJ) from Torpedo californica complexed with tacrine showed that tacrine only interacted with the CAS of AChE. However, due to its adverse effects such as acute liver toxicity and increased rates of syncope, tacrine has been gradually withdrawn from market. Since then, pharmaceutical chemistry scientists have become interested searching for AChE inhibitors able to simultaneously bind to their CAS and PAS. Several types of dual-binding-site AChE inhibitors have been developed by connecting the two interacting units through a suitable linker, which were generally derived from known AChE inhibitors either commercialized or under development13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23.

To date, donepezil (PDB entry: 4EY7) is the only dual binding site AChE inhibitor approved for the treatment of AD. The latest X-ray crystallographic structure of the complex between recombinant human AChE and donepezil reveals that the elongated structure of donepezil spans the entire length of the enzyme-active-site gorge. It has a unique orientation along the active-site gorge, extending from the CAS, at the bottom near Trp86, to the PAS at the top near Trp286. This provides a more accurate platform for further design of next-generation derivatives.

In the current work, two genistein derivatives (G1 and G2) have been discovered with strong or moderate activity against both AChE and butyrylcholinesterase (BuChE). Although G1 and G2 (Fig. 1) have quite similar structures, experimental data in this study show that the inhibitory effect of G1 against AChE was almost 80 times greater than that of G2. Since their inhibitory mechanisms against AChE are still unclear, it is of great interest to investigate why these two analogs have differing inhibitory potencies and to reveal the molecular basis for their binding to AChE. Thus, the binding mechanisms of these two inhibitors were studied by multiple approaches consisting of enzyme-kinetic analysis, molecular docking and molecular dynamics (MD) simulation26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41. This study provides a molecular basis for understanding how different configurations influence their binding affinities.

Materials and methods

In vitro inhibition studies on AChE and BuChE

AChE (E.C. 3.1.1.7) was extracted from rat cortex (Sprague Dawley). BuChE (E.C. 3.1.1.8) was obtained from human plasma (purchased from Beijing Red Cross Blood Center). 5,5′-dithiobis-(2-nitrobenzoic acid) (Ellman׳s reagent, DTNB), acetylthiocholine iodide (S-ACh), S-butyrylthiocholine chloride (BuSCh), donepezil hydrochloride and tetraisopropyl pyrophosphoramide (iso-OMPA) were purchased from Sigma-Aldrich. The AChE/BuChE activity was evaluated by detecting the hydrolysis product of the substrate following the method of Ellman. Five serial dilutions of G1, G2 and the reference compounds were utilized during the inhibition of AChE/BuChE activity.

In the AChE reaction system, the inhibition assay was carried out in 96-well plates using a Spectra Max M5 microplate reader (Molecular Devices, Sunnyvale, CA, USA). Donepezil hydrochloride was chosen as the reference compound. The reaction system consisted of 10 μL of test compounds, 30 μL of 0.05 mol/L PBS, 20 μL of AChE, 60 μL of 3.75 mmol/L of substrate S-ACh and 80 μL of 0.25 mg/mL DTNB. The assay mixture was incubated at 37 °C for 60 min and the absorbance was measured at 412 nm. Each concentration was analyzed in triplicate, and IC50 values were determined graphically from log concentration–inhibition curves.

In vitro BuChE assay was similar to the method described above. iso-OMPA served as the reference compound. The assay solution consisted of 10 μL of test compounds, 40 μL of BuChE, 70 μL of 7.5 mmol/L of substrate (BuSCh) and 80 μL of 0.25 mg/mL DTNB.

Kinetic characterization of AChE inhibition

The kinetic studies of G1 against AChE were performed using rat cortex homogenate as an AChE source. Enzyme activities were determined at 37 °C using five concentrations (50 µmol/L, 100 µmol/L, 200 µmol/L, 300 µmol/L and 500 µmol/L) of ACh in the presence or absence of three concentrations (110 nmol/L, 330 nmol/L and 1000 nmol/L) of G1. Kinetic characterization of the hydrolysis of Ach catalyzed by AChE was performed spectrometrically at 412 nm. The OD value was read for 50 min at 2-min intervals. Then the data were plotted on a Lineweaver–Burk diagram to reveal the mechanism of inhibition.

Molecular docking

The crystal structure of recombinant human AChE in complex with donepezil at a resolution of 2.35 Å was downloaded from the PDB (PDB entry: 4EY7), and was used as the initial 3D model. To obtain the docking–binding models for AChE in complexes with G1 and G2, the molecular modeling program MOE 2010 (Chemical Computing Group. Inc., Canada) was used to perform the docking process. At first, crystallographic water molecules were removed from 4EY7, while the active site pocket was defined by the ligand of donepezil. The detailed parameters were set as follows: the placement method, the first scoring function rescoring 1, and the saved poses were set to Triangle Matcher, London dG, and 30, respectively. In addition, the refinement, the second refinement scoring function rescoring 2, and the saved poses were set to forcefield, none, and 10, respectively.

After the docking parameters were set, the crystal pose of donepezil was first re-docked into the binding site pocket of AChE, and the root-mean-square deviation (RMSD) values between the docking and initial poses were calculated. Two systems were then prepared. System 1, named AChE/G1 complex, was AChE in complex with G1, whereas system 2, called AChE/G2, was AChE in complex with G2.

MD simulations

The AMBER 12 software was used to perform all the MD simulations in this study. The inhibitors were minimized using the HF/6-31G⁎ optimization in Gaussian 09, and the atom partial charges were obtained by fitting the electrostatic potentials derived by Gaussian via the RESP fitting technique in AMBER 12. The antechamber and tleap modules were used to assign the Generalized Amber Force Field (GAFF) parameters to the two systems. Ff12SB force field was assigned to the receptors. The complexes were neutralized by adding 10 sodium counterions, and were surrounded by a periodic box of TIP3P water molecules extended 8 Å between the box edges and solute surface. Simulations were performed using PMEMD.CUDA in AMBER 12. Periodic boundary conditions were applied to avoid the edge effect. The Particle Mesh Ewald (PME) method was used to calculate long-range electrostatic interactions. The hydrogen bonds were constrained using SHAKE algorithm. For the non-bonded interactions, a residue-based cutoff of 10 Å was used.

The two complexes were energy-minimized to remove possible steric stress. At first, the solvent molecules and counter ions were relaxed during a 2000-step minimization. The first 1200 steps were carried out using the steepest descent algorithm, and the remaining 800 steps were performed with a conjugate gradient algorithm. In addition, the backbone atoms of each system were fixed for 3000 cycles before 8000 cycles of fully relaxed energy minimization. The system was then slowly heated from 0 to 300 K within 50 ps. After the heating process, a further 1000 ps of equilibration at 300 K was carried out to obtain a stable density. Afterward, an unconstrained production phase was initiated and continued for 10 ns in an NPT ensemble at 1 atm and 300 K. The time step used for the MD simulations was set to 2.0 fs and the trajectory files were collected every 1 ps for the subsequent analysis. All trajectory analysis was done with the Ptraj module in AMBER 12 and examined visually using the VMD software.

MM/GBSA binding-free-energy calculations

Binding free energies were calculated using the molecular mechanics generalized-born surface area (MM/GBSA) algorithms implemented in the AMBER 12. The total binding free energy can be calculated according to the equation:where ΔEMM denotes the gas-phase interaction energy between receptor and ligand (including van der Waals energy contribution (ΔEvdw) and electrostatic energy contribution (ΔEele)); ΔGGB and ΔGSA stand for the polar and nonpolar components of the desolvation free energy, respectively; TΔS represents the conformational entropy contribution at temperature T. Here, ΔGGB was determined by the Generalized-Born approximation model, while ΔGSA was estimated based on the solvent accessible surface area (SASA) model by LCPO method: ΔGSA=0.0072×SASA.

In general, due to the relatively high computational demand of calculation of entropy in AMBER 12 and the error that may introduce significant uncertainty in the results, as well as the slight variation of substitution in genistein derivatives, the conformational entropy was not considered.

The interaction between inhibitors and each residue was computed using the MM/GBSA decomposition procedure in AMBER 12. The binding interaction of each ligand-residue pair includes three terms: ΔEvdw, ΔEele, and solvation contribution (ΔGGB+ΔGSA). ΔGGB (the polar part) was calculated using the Generalized-Born model, while ΔGSA (the non-polar part) was determined with the ICOSA method.

All energy components were computed by using the same 100 snapshots extracted from 5 to 10 ns. The free energy decomposition calculation could give insight into the mechanisms of inhibition of AChE by genistein derivatives and clarify the molecular basis for their different binding affinities.

Results and discussion

Inhibitory effect of G1 and G2 towards cholinesterase activities

As given in Table 1, donepezil hydrochloride has an IC50 value of 26 nmol/L on AChE, which is approximately equal to that of 32 nmol/L reported in the literature. G1 is a highly potent inhibitor of AChE with an IC50 value of 264 nmol/L, while G2 shows a moderate inhibitory effect toward AChE with an IC50 value of 21,210 nmol/L. Thus, the inhibitory potency of G1 is almost 80 times stronger that of G2 on AChE. To compare the selectivity of G1 and G2, an in vitro BuChE assay was also performed. The reference compound iso-OMPA exhibits an IC50 value of 2440 nmol/L with BuChE, which is comparable to that of 980 nmol/L reported in the literature. The results (Table 1) suggest both G1 (IC50=5746 nmol/L) and G2 (IC50=25,846 nmol/L) can inhibit BuChE activity moderately, while G1 (selectivity=21.8) possesses better selectivity than G2 (selectivity=1.2).

Table 1

The inhibitory effect of G1, G2 and the reference compounds on AChE and BuChE activities.

Compound	IC50 (nmol/L)a		Selectivity (AChE/BuChE)
AChE	BuChE
G1	264±24	5746±144	21.8
G2	21,210±4965	25,846±7,040	1.2
Donepezil54	26±6	–	–
ISO-OMPA55	–	2440±170	–

Note: Data are expressed as mean±SD, n=3.

Kinetic study results

Graphical analysis of the Lineweaver–Burk plot gives information about the binding mode. As shown in Fig. 2, the lines cross the third quadrant in the same point, and Vmax decreases as the concentration of G1 increases. The Lineweaver–Burk plot reveals that G1 is a typically mixed AChE inhibitor which is similar to that of donepezil55, 56. The result shows that G1 is able to interact with both CAS and PAS of AChE.

Figure 2

Steady state inhibition by G1 of AChE hydrolysis of ACh; reciprocal plots of initial velocity with substrate concentration: the plots show mixed-type inhibition for G1.

Validation of molecular docking and molecular dynamics simulation methods

To validate whether the docking software MOE was reliable for the AChE (PDB code: 4EY7) system, the crystal pose of donepezil was initially re-docked into the binding site of AChE with specific docking parameters and scoring functions. The 10 generated conformations were then compared with the crystal pose. Docking was considered acceptable if the average RMSD value of the docking poses was less than a given threshold of 2.0 Å from the crystal pose. In this study, the average RMSD value of donepezil between the docking poses and the crystal conformation was 1.53 Å. This implied that MOE could reproduce the crystal-binding model and was suitable for the AChE system. Afterward, G1 and G2 were docked into the same binding site pocket of AChE by use of the same optimum docking conditions above to obtain the starting models for subsequent MD simulations.

To explore the dynamic stability of the two AChE/inhibitor complexes and to ensure the rationality of the sampling strategy, the time-dependent RMSD values of the complex backbone atoms were calculated from the X-ray crystal structure of AChE during MD trajectories of 10 ns. The RMSD plot (Fig. 3) indicates that the conformations of the AChE/G1 complex achieve equilibrium around 3200 ps and fluctuate around 1.3 Å. However, for the AChE/G2 complex, the equilibrium time is around 5000 ps and the conformations fluctuate around 1.6 Å. Both trajectories are stable after 5000 ps, so it is reasonable to do the binding free energy calculation and free energy decomposition based on the snapshots extracted from 5 to 10 ns.

Figure 3

Root-mean-square deviation (RMSD) of the backbone atoms (CA, N, C) of the AChE/G1 and AChE/G2 complexes.

For better understanding on how G1 and G2 influence the binding mode with AChE, the structural changes of two complexes were examined. To visualize the result clearly, 10 structures (snapshots) were overlaid during the MD simulations, of which coordinates were saved after every 1 ns from 1 to 10 ns. As shown in Fig. 4a, the 10 conformations of G1 undergo little movement, which means that G1 is stable throughout the MD simulation. On the other hand, in Fig. 4b, the 10 conformations of G2 undergo a large amount of movement during the MD simulations. Therefore, this analysis suggests that AChE/G1 complex is more stable than AChE/G2 complex.

Figure 4

The superposition of 10 snapshots from 1 to 10 ns for AChE/G1 complex (a) and AChE/G2 complex (b).

Binding free energies calculation

In MM/GBSA calculations, the affinity of one inhibitor binding to the receptor could be estimated using the snapshots from a trajectory of the complex. The ΔGexp was estimated by the equation ΔGexp≈RT ln(IC50). These values, with the values of absolute binding free energies of G1 and G2, obtained using the MM/GBSA technique, are given in Table 2. The descending trend in ΔGpred for AChE/G2 and AChE/G1 complexes is consistent with that in ΔGexp.

Table 2

The components of the binding free energy (kcal/mol) for the AChE complexes with G1 or G2, determined by using the MM/GBSA method.

Inhibitor	ΔEvdw	ΔEele	ΔGGB	ΔGSA	ΔGpred	ΔGexp
G1	−54.8±2.8	−11.6±2.2	29.1±2.2	−6.6±0.2	−43.9±2.9	−9.0
G2	−54.8±2.4	−3.9±3.2	26.0±2.8	−6.2±0.2	−38.9±2.6	−6.4

To get a better view of which energy term has more impact on the binding affinity, the four individual energy components (ΔEvdw, ΔEele, ΔGGB and ΔGSA) were carefully compared. From Table 2, it is evident that the ΔEvdw is dominant in ΔGpred. The ΔEvdw of AChE/G1 (−54.8 kcal/mol) complex is equal to AChE/G2 (−54.8 kcal/mol). Obviously, the ΔEvdw part cannot account for the activity difference between G1 and G2. However, the ΔEele of AChE/G1 complex (−11.6 kcal/mol) is stronger than that of AChE/G2 (−3.9 kcal/mol). Considering the negative effect generated by ΔGGB, the net electrostatic contributions (ΔEele+ΔGGB) of AChE/G1 and AChE/G2 complexes are 17.5 and 22.1 kcal/mol, respectively, which is unfavorable for binding in the two complexes. Furthermore, the difference in the net electrostatic contribution (ΔEele+ΔGGB) between G1 and G2 is −4.6 kcal/mol. In addition, both of ΔGSA for the two complexes (−6.6 kcal/mol and −6.2 kcal/mol) are approximately equivalent. Therefore, the net electrostatic contributions (ΔEele+ΔGGB) play the most important role in differentiating the activity between G1 and G2.

Residue-based-energy decomposition using the MM/GBSA method

To further delineate the detailed mechanism of the AChE/inhibitor interactions, residue-based-energy decomposition was performed to evaluate the effects of energy on the contributions of each residue in the binding site pocket by using the MM/GBSA method. Fig. 5 shows the energy decomposition values for key residues in the two complexes.

Figure 5

The inhibitor–residue interaction spectrum for AChE/G1 and AChE/G2 complexes.

Generally, if the interaction energy between a residue and a ligand is lower than −0.8 kcal/mol, the residue is regarded as an important residue in the molecular recognition of that ligand. In the AChE/G1 complex, the major favorable energy contributions (−3.02 to −0.94 kcal/mol) originate predominately from Tyr124 (−3.02), Trp341 (−2.18), Trp286 (−2.08), Glu292 (−1.49), Val294 (−1.28), Phe338 (−1.04) and Tyr72 (−0.94). In the AChE/G2 complex, the residues which give the most negative binding free energies (kcal/mol) are Trp286 (−2.79), Tyr341 (−2.43), Thr83 (−1.23), Gln291 (−0.93), Trp86 (−0.92), Tyr124 (−0.90) and Glu292 (−0.88). As shown in Fig. 5, the two inhibitors have similar interaction patterns, which means strong or moderate interactions with residue Tyr124, Trp286, Glu292 and Try341. Comparison between G1 and G2 shows that the interactions of G1 with Tyr124, Glu292, Val294 and Phe338 (−3.02, −1.49, −1.28 and −1.04 kcal/mol) are stronger than those of G2 (−0.90, −0.88, −0.20 and −0.36 kcal/mol). Therefore, it is conceivable that the difference between the inhibitory effect of G1 and G2 originates from these key residues (Tyr124, Glu292, Val294 and Phe338).

Since ΔEvdw is dominant in the total binding free energy (ΔGpred), it is of interest to determine if van der Waals interactions determine the different biological activities between G1 and G2. As shown in Supplementary Fig. 1, the differences in the van der Waals interactions between G1 and G2 were reported for the amino acids Trp286, Val294, Phe338 and Tyr341. The van der Waals interactions of G1 with Val294 and Phe338 (−0.97 and −1.93 kcal/mol) are stronger than those of G2 (−0.09 and −0.64 kcal/mol). However, the interactions of G1 with Trp286 and Tyr341 are weaker (−3.48 and −2.81 kcal/mol) than those of G2 (−4.28 and −3.68 kcal/mol), which compensates for the differences made by van der Waals interactions.

As presented in Supplementary Fig. 2, electrostatic interactions (ΔEele) between the inhibitors and the important residues were investigated to explain the different inhibitory activities between G1 and G2. It is obvious that the ΔEele of G1 with Tyr124, Glu292, Phe338 and Tyr341 (−1.97, 0.06, −0.21 and −0.44 kcal/mol) is much stronger than those of G2 (−0.94, 0.66, 0.39 and 0.33 kcal/mol). Therefore, the differences obtained from electrostatic interactions can explain the significant difference between G1 and G2.

As shown in Supplementary Fig. 3, the polar contributions of desolvation (ΔGGB) were also determined. The interaction is almost opposite to the electrostatic interaction. On the whole, the net electrostatic interaction (ΔEele+ΔGGB) opposes binding, which is consistent with the discussion in the section on binding free energies calculation.

Binding mode of the AChE/G1 and AChE/G2 complexes

Results above proposed an explanation as to why G1 was more active than G2 in terms of free energy. In fact, the energy was decided by the structure. Therefore, we analyzed whether slight structural modifications on G1 could influence the binding and eventually change the interaction with AChE. Here the average structure was based on the 100 snapshots (from 5 to 10 ns) obtained, and the binding modes in the active site pocket of AChE with G1 and G2 are displayed in Fig. 6.

Figure 6

The binding modes of AChE/G1 (a) and AChE/G2 (b) complexes.

As shown in Fig. 6a, G1 extends deeply into the two main binding sites of AChE. The CAS and PAS are connected by a deep, hydrophobic gorge. The piperidine ring of G1 enables hydrophobic/van der Waals interactions with Trp86 and Phe338 in the CAS, while the flavone ring of G1 stacks against Trp286 and Tyr341 in the PAS; thus G1 is sandwiched between the rings of Trp286 and Tyr341. Meanwhile, the oxygen atom at C7 in flavone ring can form a strong hydrogen bond of 2.8 Å with the –OH of Tyr124, which is consistent with decomposition analysis of the electrostatic interaction. Moreover, G1 can also make hydrophobic/van der Waals interactions with Tyr72, Glu292 and Val294.

As shown in Fig. 6b, G2 has a similar binding pattern to G1 except for some key differences. The flavone ring of G2 can also stack against Trp286 and Tyr341 in the PAS and make a hydrophobic/van der Waals interaction with Trp86 in the CAS. Because of the long distance (3.5 Å) between the oxygen atom at C7 and the –OH of Tyr124, they could not form a hydrogen-bonding interaction. In order to clearly observe the differences between the two binding modes, the two complexes were superimposed. As presented in Fig. 7, compared with G1, the methyl group substitution at C16 induces a large movement of the piperidine ring of G2 as demonstrated by MD. Therefore, some hydrophobic/van der Waals interactions (such as Tyr72, Val294 and Phe338) have disappeared.

Figure 7

The superposition of AChE/G1 (yellow) complex over AChE/G2 (cyan) counterparts. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Conclusions

Two genistein derivatives (G1 and G2) have been discovered as inhibitors of AChE and BuChE in this investigation. The absence (G1) or presence (G2) of a methyl group resulted in a distinct difference in their inhibitory activity toward AChE. The present study provided insights into the two inhibitor-bound structures of AChE. We analyzed the binding mechanisms of both inhibitors to AChE using multiple approaches including enzyme-kinetics analysis, molecular docking and MD simulation.

The results obtained by the kinetic study revealed that G1 was able to interact with both the CAS and PAS of AChE. MM/GBSA calculations predicted binding free energies of −43.9 kcal/mol (G1) and −38.9 kcal/mol (G2), which were consistent with our experimental results. It was observed that ΔEele, ΔEvdw and ΔGSA were the main driving forces for inhibitor binding, whereas the net electrostatic contribution (ΔEele+ΔGGB) difference was the leading reason for distinguishing the binding affinities of these two inhibitors.

The structure analysis revealed that the flavone rings of the two inhibitors could stack against Trp286 and Tyr341 in the PAS, and their piperidine rings could enable hydrophobic/van der Waals interaction with Trp86 in the CAS. This may explain why genistein derivatives with a piperidine ring exhibit inhibition of AChE. In the binding mode of G1, van der Waals and electrostatic/H-bond interactions are the driving forces for AChE recognition. Unlike G1, the oxygen atom at C7 in G2 could not form a hydrogen bond with the –OH of Tyr124, and some hydrophobic/van der Waals׳ interactions (such as Tyr72, Glu292, Val294 and Phe338) had disappeared, making the inhibitory potency of G1 80 times higher than that of G2.

In summary, all the results indicated that enzyme-kinetic analysis, molecular docking and MD simulation could provide an alternative way to reveal the molecular mechanism of inhibition of AChE by genistein derivatives, and thus contribute to a better understanding of the therapeutic potential of the two inhibitors.