Molecular Modeling Studies of Piperidine Derivatives as New Acetylcholinesterase Inhibitors against Neurodegenerative Diseases
Neurodegenerative disorders are related to the progressive loss of structure or function and, eventually, death of neurons. These processes are responsible for diseases like Parkinson’s, Alzheimer’s, and Huntington’s, and the main molecular target for the drug design against these illnesses today is the enzyme acetylcholinesterase (AChE). Following this line, in the present work, we applied docking techniques to study some piperidine derivative inhibitors of AChE and further propose structures of six new AChE inhibitors as potential new drugs against neurodegenerative disorders. The best inhibitor proposed was submitted to additional molecular dynamics simulations steps.
Skeletal muscle ontogenesis is a complex and multistep process during which myoblasts, the skeletal muscle progenitors, proliferate, migrate, and fuse into myotubes . The latter are then innervated, and the nerve endings govern the final phases of maturation of myotubes into the functional skeletal muscle fibres responsive to nerve inputs . Regular transmission of these inputs is directly related to the action of acetylcholinesterase (AChE), an enzyme present in the neuromuscular junctions and responsible for the hydrolysis of the neurotransmitter acetylcholine. In some neurodegenerative disorders, like myastimiahenia gravis or Parkinson’s disease, inhibition of AChE is a fundamental therapeutic step [1, 3].
Several clinical studies on AChE inhibitors, such as physostigmine and tacrine, have revealed interesting results. However, side effects and oral activity restrict the application of these agents as effective drugs. In this context, recent findings show that piperidine derivatives are promising AChE inhibitors able to overcome the unfavorable side effect profile and poor pharmacokinetics related to the aforementioned compounds. Despite great importance, the study of the interaction between piperidine derivatives and AChE has received remarkably little attention. The interactions of AChE with its inhibitors, however, are not fully understood yet and still deserve investigation. The drug design of rational therapy requires insight into the molecular mechanisms underlying the AChE-inhibitor interactions under physiological conditions in order to propose structures of new and more efficient inhibitors. In this line, molecular mechanics based methods involving docking studies and also molecular dynamics simulations are suitable tools to adjust ligands at target sites and to estimate interaction energy (affinity) . Nowadays, that is a well-established technique applied to numerous cases [5, 6]. Thus, the aim of this work is to investigate the interactions among the potential inhibitors and AChE by docking and molecular modeling techniques in order to contribute to the elucidation of their action mechanisms.
2.1. Ligands Data Set
The compounds for analysis (Table 1) were selected from a series of piperidine derivatives, known as potent inhibitors of AChE, with (half maximal inhibitory concentration) values in the μM range. The 3D structures of each compound in the neutral forms were constructed using the Spartan software [7, 8]. Each structure was geometry-optimized without any restriction in vacuum, and the partial atomic charges were assigned using the AM1 semiempirical method [9–11].
2.2. Docking Studies
The 3D structure of AChE from Torpedo californica including waters of crystallization and donepezil was obtained from the Protein Data Bank (PDB—code 1EVE) . The compounds in Table 1 were docked into the AChE binding site using the Molegro Virtual Docker (MVD) , a program for predicting the most likely conformation of how a ligand will bind to a macromolecule, and some relevant residues of the protein were considered flexible during the docking simulation. The MolDock scoring function used by MVD is derived from the piecewise linear potential (PLP), a simplified potential whose parameters are fit to protein-ligand structures and binding data scoring functions and further extended in generic evolutionary method for molecular DOCK with a new hydrogen bonding term and new charge schemes. The docking scoring function, , is defined by the following energy terms: where is the ligand-protein interaction energy: The term is a “piecewise linear potential” using two different sets of parameters: one set for approximating the steric (van der Waals) term between atoms and another stronger potential for hydrogen bonds. The second term describes the electrostatic interactions between charged atoms. It is a Coulomb potential with a distance-dependent dielectric constant given by . The numerical value of 332.0 adjusts the units of the electrostatic energy to kilocalories per mole ; is the internal energy of the ligand:
The first term (double summation) is between all atom pairs in the ligand, excluding atom pairs which are connected by two bonds. The second term is a torsional energy term, where is the torsional angle of the bond. The average of the torsional energy bond contribution is used if several torsions could be determined. The last term, , assigns a penalty of 1000 if the distance between two heavy atoms (more than two bonds apart) is less than 2.0 Å, ignoring unfeasible ligand conformations . In summary, these functions are used to automatically superimpose a flexible molecule onto a rigid template molecule. The docking search algorithm used in MVD is based on an evolutionary algorithm, where interactive optimization techniques inspired by Darwinian evolution theory and a new hybrid search algorithm called guided differential evolution are employed. The guided differential evolution algorithm combines the differential evolution optimization technique with a cavity prediction algorithm during the search process, which allows for a fast and accurate identification of potential binding modes (poses). The active site exploited in docking studies was defined, with a subset region of 10.0 Å from the center of the ligand. The interaction modes of the ligand with the enzyme active site were determined as the highest energy scored protein-ligand complex used during docking, and the conformers of each compound were mostly associated with bioactive conformations of donepezil [14–16].
2.3. QM/MM Calculations
In order to investigate the influence of electronic effects on the relative binding energy, we have performed QM calculations of piperidine derivatives in the AChE active site. The model system used was the ligand and some residues of amino acids from the active site. The results were performed by calculating the relative binding energy at the QM/MM level from docking structures. In all cases, the QM calculations were acquired at the density functional theory (DFT) level with the Gaussian 03 code . DFT methods have shown an excellent performance for medium and large systems and have also proven to be appropriate for biomolecules . Calculations were performed using the generalized gradient approximation functional proposed by Gustin et al. . The optimized structures were improved by including single-point energies using the 6–31 G (d, p) basis set, zero-point energy, and thermal corrections (at 298.15 K and 1 atm) estimated at the B3LYP/6-31 G (d) level.
2.4. Molecular Dynamic Simulations
The best pose of compound 2 inside AChE, obtained from the docking studies, was submitted to additional MD simulation steps carried out using the GROMACS 4.5.4 package [20, 21] with the force field OPLS/AA [22, 23]. To obtain the parameters and topologies for the referred compound, the AnteChamber PYthon Parcer InterfacE (ACPYPE)  was used. The software MKTOP [20, 21] was used to generate the parameters that were not given by ACPYPE . OPLS/AA [22, 23] force field parameters for the ligands were used, and the atomic partial charges were calculated by semiempirical quantum chemistry SQM [22, 23] (via Antechamber) according to AM1-BCC parameters (parameterized to reproduce HF/6–31 G RESP charges). SQM was modified to include six decimals of precision instead of the default 3.
The dynamic behavior of the complex AChE-compound 2 was simulated in a cubic box of approximately 400 nm3, containing around 8000 water molecules of the SPC type, using periodic boundary conditions [24, 25]. In this study, for both pressure and temperature, we used Berendsen’s weak coupling algorithm scheme which maintained temperature and pressure at 300 K and 1 bar. This system was submitted to four energy minimization steps as follows: steepest descent with restrained position of the ligands and convergence criterion of 50.00 kcal mol−1 nm−1, followed by steepest descent without restriction, conjugate gradients, and, finally, quasi-Newton Raphson until an energy of 10.00 kcal mol−1 nm−1.
The minimized complex was then submitted to MD simulations in two steps. Initially, we performed 500 ps of MD at 300 K with PR for the entire system, except the water molecules, in order to ensure a balance of the solvent molecules around protein and ligand. Subsequently, 20 ns of MD was performed without any restriction, using 2 fs of integration time and a cutoff of 10 Å for long-distance interactions. The electrostatic interactions were calculated using the particle mesh Ewald method (PME), and a total of 1,000 conformations were obtained during each simulation. In this step, the lists of pairs (pairlists) were updated every 500 steps; all Arg and Lys residues were assigned with positive charges, and the residues Glu and Asp were assigned with negative charges.
To analyze the structures generated after the optimization and MD steps, we used the Visual Molecular Dynamics (VMD)  and Swiss-Pdbviewer  programs. Plots of total energy variation, distance, random mean square deviation (RMSD), and H-bonds formed during the MD simulation were generated with the Origin program . Qualitative spatial RMSD pictures were generated in MolMol program , and the figure of the MD simulation frames was generated in the PyMOL program .
3. Results and Discussion
3.1. Redocking Procedure
The docked binding mode was used to establish a link between the MolDock scoring function, structural properties of piperidine analogs, and their biological activity against the AChE.
Initially, we have performed a redocking calculation. The ligand binding process was validated by “redocking” of the cocrystallized donepezil with Torpedo californica AChE. In docking simulations, we expect that the best results generate RMSD values below 2.0 Å, when compared to crystallographic structures . We verified that donepezil became bound in the AChE crystal structure (Figure 1), similarly to preexisting cocrystallized ligand (RMSD = 0.62).
3.2. Docking Studies
The potential binding sites of AChE were calculated using Molegro cavity prediction algorithm. A cavity of 230.40 A3 (surface = 574.72 A2) was observed in the AChE and the side chains of amino acid residues from 6 Å of the cavity surface were considered flexible.
In our docking studies, the following parameters were calculated: (a) hydrogen bonding energy between the ligand and the protein (); (b) internal energy values of ligand (); (c) short-range ( Å) electrostatic protein-ligand interaction energy (); (d) long-range ( Å) electrostatic protein-ligand interaction energy (); and (e) steric interaction energy between protein and ligand (). Multiple linear regression (MLR) was used to model the relationship between the independent energy term values and values by fitting a linear equation to the observed data. The predicted was obtained by the following equation: where , , , , , , and are coefficients obtained from the multivariable regression.
The energy term values of the best docked conformations correlated with inhibition activity showed a coefficient of correlation () equal to 0.83 (see (5)). An greater than 0.7 indicates a high correlation between the data [15, 16]. This suggests that the binding conformations and binding modes of the AChE inhibitors are satisfactory. Equation (5) was established with 8 compounds (donepezil and compounds 1–7):
It should be kept in mind that lower negative energy values (most negative) represent a higher stability of complexes between ligand and protein. has great weight in the potency of the compounds since the magnitude of its regression coefficient is 0.042 (see (5)). Thus, the negative coefficients correspond to favorable interactions between the ligand and the amino acid residues in the protein active site. In other words, these interactions increase the potency of the inhibitors. For example, donepezil is the most potent compound and is shown to be the most stable (−44.13 kcal·mol−1) (Table 2). This is in good agreement with experimental data.
Equation (5) was used to determine the activity of compounds 9–14 (Table 1). Compounds 1, 2, 3, 4, 7, and 8 present the same substituent in relation to compounds 9, 10, 11, 12, 13, and 14, respectively. The difference is in the benzyl/benzhydryl ring. Table 3 compares values of compounds 1–14. Generally all compounds with benzyl group were less potent than compounds with benzhydryl group.
By analyzing the interaction energies from the docking calculation (Table 2) as well as the hydrogen bond formed between the compounds and AChE, it was observed that (i) compounds 1, 3, and 8 interact with Phe288; (ii) compounds 2 and 7 interact with Tyr121; (iii) compound 1 interacts with Tyr70. This fact could justify the favorable of −5.91 kcal·mol−1. Cation-π and π-staking interactions were observed between the benzhydryl ring and Trp84, Glu199, Ser200, and/or His440. compound 2 exhibited the lowest (more stable complexes). Figure 2 shows compound 2 docked into the AChE active site.
By correlating the docking results obtained with experimental data, it is possible to suggest the proposal of the structures of seventeen derivatives, shown in Table 4, as potential AChE inhibitors. These compounds were further evaluated using docking simulations at the AChE active site. The predicted values given in Table 4 for the proposed compounds suggest that just compound G, in this set of compounds, is more potent than donepezil. The theoretical inhibitory potencies obtained for compounds D, E, J, K, and N were, however, higher than for other experimentally more active compounds. Figure 3 shows compound G docked into the AChE active site. It interacts with Phe288 through the hydrogen bond. Trp84 forms a parallel aromatic stacking interaction with compound G, while Phe330 forms edge-to-face aromatic stacking contacts.
3.3. QM/MM Calculations
Now, in order to evaluate the electronic effects on the relative binding energy between ligand and protein, we performed QM/MM calculations. The QM region for AChE was optimized with the two-layer ONIOM method , implemented in the Gaussian 03 program . In these QM/MM calculations, a specified region around the active center was calculated with a QM method, while the rest of the protein was treated at an MM level. At the QM/MM border, atoms in the MM region, bound to an atom in the QM region, were replaced by hydrogen atoms at the QM-level part of the QM/MM calculation. The QM region includes the ligand and its surrounding amino acid residues. We decided to obtain the QM region by surrounding the ligands with a sphere of 15 Å centered, including some amino acid residues from the active site . The QM calculations were performed with the B3LYP method, which consists of the Slater exchange, the Hartree-Fock exchange, the exchange functional of Becke, the correlation functional of Lee, Yang, and Parr (LYP), and the correlation functional of Vosko et al. . The method selected for MM calculations was the Amber force field (Amber96). From Table 2 (single-point calculations with QM/MM method), the most stable orientations for the compound were used. The relative binding energies at the QM/MM level showed a very good agreement with the docking energy calculation. This means that the orientations and hydrogen bonds of piperidine derivatives into the AChE active site can be evaluated by docking calculations using molecular mechanics methods.
3.4. MD Simulation
The best pose of compound 2 obtained in the docking studies was submitted to 20 ns of MD simulations inside AChE in order to corroborate the docking results and assess the dynamical behavior of this compound. The energy plot (data not shown) for the system simulated showed that the total energy tends to stability since the first ps of simulation.
In the plot of temporal RMSD, presented in Figure 4, it is also possible to observe the equilibration of the system since the first ps and that the RMSD values never surpassed 0.40 nm during the simulation. The spatial RMSD for the system is shown in Figure 5. As can be seen, the regions with major fluctuations during the simulation (larger tubes in the picture) correspond to the extremities and loops. On the other hand, the residues in the active site region and in the α-helix and β-sheets presented lower fluctuations (lower thickness in the tubes). The results above attest to the stability of the system during the MD simulation.
Regarding the dynamical behavior of compound 2 during the simulation, it was observed that this ligand remained stable inside AChE surrounded by residues Tyr70, Trp84, Tyr121, Glu199, Ser200, Phe288, and His440 as can be seen in the superposition of frames in Figure 6. No H-bond was observed between compound 2 and AChE during the simulation time and all interactions observed were hydrophobic, involving mainly Trp84, Glu199, Ser200, and His440 residues. These observations corroborate with the docking results. A movie showing the dynamical behavior of compound 2 inside AChE can be accessed as Supplementary Material (see movie in Supplementary Material available online at http://dx.doi.org/10.1155/2013/278742) of this paper.
In the present work, seventeen new AChE inhibitors were proposed as potential new drugs against neurodegenerative disorders. Our results suggest that compound G had presented better values than donepezil. compounds D, E, J, K, and N showed theoretical inhibitory potencies better than experimentally tested drugs. We believe that these results validate the proposals as suitable leads to the drug design against neurodegenerative disorders. Experimental studies should be performed on these compounds in order to make room for further refinements on their structures that will eventually lead to new and more effective AChE inhibitors.
The authors thank the Brazilian agencies FAPEMIG, CAPES, and CNPq for funding part of this work.
Movie, created with the program VMD, showing the dynamical behavior of compound 2 inside AChE during the 20 ns of simulation. For better clarity only residues Trp84, Glu199, Ser200, and His440 (involved in hydrophobic interactions with compound 2) are shown. No H-bond interaction was observed during the simulation time.
K. Shimizu, H. Fujita, and E. Nagamori, “Evaluation systems of generated forces of skeletal muscle cell-based bio-actuators,” Journal of Bioscience and Bioengineering, vol. 115, pp. 115–121, 2013.View at: Google Scholar
K. S. Matos, E. F. F. da Cunha, A. da Silva Gonçalves et al., “First principles calculations of thermodynamics and kinetic parameters and molecular dynamics simulations of acetylcholinesterase reactivators: can mouse data provide new insights into humans?” Journal of Biomolecular Structure and Dynamics, vol. 30, pp. 546–558, 2012.View at: Publisher Site | Google Scholar
M. S. Caetano, T. C. Ramalho, T. G. Vieira, A. D. Goncalves, D. T. Mancini, and E. F. F. da Cunha, “Bonding, structure, and stability of clusters: some surprising results from an experimental and theoretical investigation in gas phase,” Journal of Chemistry, vol. 2013, Article ID 362894, 8 pages, 2013.View at: Publisher Site | Google Scholar
SPARTAN, “Version 10.0,” Wavefunction, Irvine, Calif, USA.View at: Google Scholar
D. C. Young, Computational Chemistry: A Practical Guide For Applying Techniques To Real-World Problems., Wiley-Interscience, New York, NY, USA.
C. J. Cramer and D. G. Truhlar, “General parameterized SCF model for free energies of solvation in aqueous solution,” Journal of the American Chemical Society, vol. 113, no. 22, pp. 8305–8311, 1991.View at: Google Scholar
C. J. Cramer and D. G. Truhlar, “An SCF solvation model for the hydrophobie effect and absolute free energies of aqueous solvation,” Science, vol. 256, no. 5054, pp. 213–217, 1992.View at: Google Scholar
M. Rooman, E. Cauët, J. Liévin, and R. Wintjens, “Conformations consistent with charge migration observed in DNA and RNA X-ray structures,” Journal of Biomolecular Structure and Dynamics, vol. 28, no. 6, pp. 949–954, 2011.View at: Google Scholar
M. J. Frisch, G. W. Trucks, H. B. Schlegel et al., Gaussian, Pittsburgh, Pa, USA, 1998.
A. A. S. T. Ribeiro, B. A. C. Horta, and R. B. De Alencastro, “MKTOP: a program for automatic construction of molecular topologies,” Journal of the Brazilian Chemical Society, vol. 19, no. 7, pp. 1433–1435, 2008.View at: Google Scholar
C. Oostenbrink, A. Villa, A. E. Mark, and W. F. van Gunsteren, “A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force-field parameter sets 53A5 and 53A6,” Journal of Computational Chemistry, vol. 25, no. 13, pp. 1656–1676, 2004.View at: Publisher Site | Google Scholar
D. van der SPOEL, A. R. Buuren, and E. APOL, GROMACS User Manual Version 3. 0. Groningen, University of Groningen, Department of Biophysical Chemistry, 2001.
K. S. Matos, D. T. Mancini, E. F. F. da Cunha, K. Kuca, T. C. C. Franca, and T. C. Ramalho, “Molecular aspects of the reactivation process of acetylcholinesterase inhibited by cyclosarin,” Journal of Brazilian Chemical Society, vol. 22, pp. 1999–2004, 2011.View at: Google Scholar
H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, and J. Hermans, Interaction Models for Water in Relation to Protein Hydration, Reidel, Dordrecht, The Netherlands, 1981.
W. L. DeLano and S. Bromberg, PyMOL USer’S Guide, DeLano Scientific LLC, San Francisco, Ca, USA, 2004.
T. C. Ramalho, E. F. F. da Cunha, and R. B. De Alencastro, “Solvent effects on13C and15N shielding tensors of nitroimidazoles in the condensed phase: a sequential molecular dynamics/quantum mechanics study,” Journal of Physics Condensed Matter, vol. 16, no. 34, pp. 6159–6170, 2004.View at: Publisher Site | Google Scholar
S. H. Vosko, L. Wilk, and M. Nusair, “Accurate spin-dependent electron liquid correlation energies for local spin-density calculations—critical analysis,” Canadian Journal of Chemistry, vol. 58, pp. 1200–1211, 1980.View at: Google Scholar