Molecular dynamics of the interaction of pralidoxime and deazapralidoxime with acetylcholinesterase inhibited by the neurotoxic agent tabun

ARTICLE

Molecular dynamics of the interaction of pralidoxime and deazapralidoxime with acetylcholinesterase inhibited by the neurotoxic agent tabun

Arlan da S. Gonçalves^I; Tanos C. C. França^I; Alan Wilter^II; José D. Figueroa-Villar^{*, I}

^IDepartamento de Química, Instituto Militar de Engenharia, Pça General Tibúrcio 80, 22290-270 Rio de Janeiro-RJ, Brazil

^IIBiblioteca - Memória Técnico-Científica, Instituto Nacional de Pesquisas Espaciais (INPE/MCT), Av. dos Astronautas, 1758 - Jardim da Granja, 12227-010 São José dos Campos-SP, Brazil

ABSTRACT

Efficient acetylcholinesterase reactivators are fundamental for the development of antidotes against poisoning by neurotoxic pesticides and chemical warfare agents. However, the mechanism of the reactivation reaction and the structural characteristics of the known reactivators are poorly understood. In order to study the dynamic behavior and the effect of the antidote net charge in the reactivation of this enzyme, we carried out a molecular dynamics study of human acetylcholinesterase inhibited by tabun in complex with the antidote pralidoxime and with its deaza analogues in the neutral and anionic forms. Results show that the positive charge of pralidoxime is important for its admission and permanence inside the active site. Also, the analogues, unlike pralidoxime, when forced inside the active site, move away from the phosphorilated serine residue of the enzyme and are repelled by the electrostatic potential at the entrance of the channel that conducts to the active site.

Keywords: acetylcholinesterase, molecular dynamics, tabun, antidotes, neurotoxic agents

RESUMO

Reativadores eficientes de Aceticolinesterase são fundamentais para o desenvolvimento de antídotos contra o envenenamento por pesticidas neurotóxicos e agentes de guerra química. Todavia, o mecanismo da reação de reativação e as características estruturais dos reativadores conhecidos são pouco compreendidos. Com o objetivo de estudar o comportamento dinâmico e o efeito da carga líquida do antídoto na reativação desta enzima, foi conduzido um estudo por dinâmica molecular da acetilcolinesterase humana inibida por tabun em complexo com o antídoto pralidoxima e com seu análogo deazapralidoxima nas formas neutra e aniônica. Os resultados mostraram que a carga positiva da pralidoxima é importrante para sua admissão e permanência dentro do sítio ativo. Além disso, os análogos, diferente da pralidoxima, quando colocados dentro do sítio ativo, se distanciam do resíduo serina fosforilado da enzima e são repelidos pelo potencial eletrostático na entrada do canal que conduz ao sítio ativo.

Introduction

The intensive use of neurotoxic organophosphorous compounds as pesticides in agriculture, as well as their potential use as mass destruction agents in chemical warfare, has attracted attention to the development of efficient antidotes for this type of poisoning.^1,2 However, the knowledge on the appropriate treatment for patients exposed to this kind of compounds is limited to few groups of physicians around the world.^1,2

One particularly important family of lethal tactical warfare chemicals is the group known as the nerve agents, which are closely related in chemical structure and biological action to many commonly used organophosphorous insecticides, but which are much more lethal.^1,2

The nerve agents are esters of phosphoric acid and are potent inhibitors of acetylcholinesterase, a fundamental enzyme for ending nervous impulses. These compounds inhibit all acetylcholinesterases, including the human enzyme (HuAChE), by phosphorylating a serine hydroxyl group (Ser203 in HuAChE), which is directly responsible for the hydrolysis of the neurotransmitter acetylcholine. This reaction occurs very rapidly and can lead to irreversible inhibition by a process called aging.³ Before aging, the inhibition of AChE can be reversed by dephosphorilation of the serine residue by a nucleophile, usually an oxime. In fact, it is believed that the hydroxyl group of the oxime carries out a nucleophilic attack on the phosphorylated serine residue, removing the phosphate moiety and reactivating the enzyme,³ as shown in Figure 1. Accordingly, the standard treatment for intoxication with nerve agents involves the administration of an anticholinergic drug, usually intravenous atropine, and a cationic oxime, like pralidoxime (2-PAM).³

In a previous work⁴ it was shown, by ab initio calculations, that considering the two forms of 2-PAM that are available in physiologic conditions, the protonated form is the most flexible and displays a greater negative charge at the oxime oxygen. Similar results were obtained for the bicationic oxime HI-6 (1-(2-hydroxyimino-1-methylpyridine)-1-(4-carboxyiminopyridine)-dimethyl ether hydrochloride).⁵

In this work, we have used docking and molecular dynamics studies to determine the role of the net charge of oximes as antidotes for tabun (GA) inhibition in HuAChE. For this, we carried out molecular dynamics simulations for 2-PAM and its deaza analogue in its neutral form (DZP) and its anionic form (DZP_anion), inside and outside the active site of HuAChE inhibited with GA.

Methodology

The crystallographic structure of HuAChE refined at 2.76 Å of resolution and complexed with fasciculin-II (PDB code 1B41) obtained by Kryger et al.,⁶ in early 1999, using X-ray Diffraction, left some missing residues between the amino acids Pro492 and Pro495 and Pro258 and Asn265. These missing residues were modeled using the complete FASTA format sequence of HuAChE fitted over the incomplete 3D structure without the Fasciculin-II, by means of the Swiss-PdbViewer program.⁷ Further, this alignment was submitted to the Optimize mode of the Swiss-Model server^8,9 that generated the initial complete model of the target enzyme. This model was minimized using the GROMOS 87 force field,¹⁰ implemented in the computational package GROMACS 3.2.1^10,11 until a gradient of 418 kJ mol^-1 nm^-1 (10 kcal mol^-1 Å^-1). The steepest descent algorithm was employed during this minimization procedure and the minimized model was further submitted to the BIOTECH server¹² for validation.

2-PAM, DZP and DZP_anion were manually fitted at a conformation proposed by Castro et al.⁴ as the best one for a good interaction with HuAChE inhibited with GA. After that, the coordinates of neutral DZP, DZP_anion, 2-PAM and SGA (Serine 203 bonded with GA) were submitted to the Dundee PRODRG Server¹³ in order to obtain the topologies needed for posterior MD simulations with the GROMACS 3.2.1 package.^10,11 As the GROMOS 87 force field¹⁰ does not present topologies for other moieties than aminoacids or nucleotides, the parameters for the SGA atoms, bonded lengths, angles and dihedral angles were added to their databank (ffgmx.rtp file)¹⁴ in order to create a residue having a covalent bond between GA and Ser203 (aminoacids.dat file).¹⁴ Charge calculations for 2-PAM, DZP, DZP_anion and SGA atoms were carried out using the software GAMESS US¹⁵ at the Hartree-Fock/6-31G^* level and using the CHELPG algorithm. These values are shown in Table 1. After parametrization and charge calculations, the HuAChE molecules with Ser203 replaced by SGA were manually docked with 2-PAM, DZP and DZP_anion inside and outside the active site, giving origin to nine enzyme-tabun-antidote systems, which were called HuAChE-GA-ANT systems. All these systems were further submitted to energy minimizations and MD simulations steps. The structures of GA, 2-PAM, DZP and DZP_anion are presented in Figure 2.

Previous to the MD simulations each system was placed inside a 722.86 nm³ cubic box containing 21,332 water molecules and their energy minimized by the steepest descent algorithm, until reaching a gradient of 418 kJ mol^-1 nm^-1 (10 kcal mol^-1 Å^-1). The subsequent MD steps were carried out according to the following procedure: first, 50 ps of molecular dynamics at 300 K, with all atoms having positions restrained except for the water molecules inside the box. This initial MD was necessary in order to allow for the equilibration of the solvent molecules around the protein. Then, there were carried out the full MD simulations of 1000 ps at 300K with no restrictions, using 2 fs of integration time, PME^16,17 for long-range electrostatic interactions, with 1.4 nm in the real space, and a cut-off of 1.4 nm for long-range van der Waals interactions. As a whole, 50 conformations were stored during each simulation. The pair-lists were updated every 1000 steps. For each compound there were carried out three MD simulations. The first one with an initial position of the compound inside the inhibited enzyme active site, at 0.3 nm between the phosphorous atom of SGA and the oxygen atom of the oximes, the second also inside the active site but with an O-P distance of 0.5 nm, and the last one with the compound at the entrance of the well that leads to the active site (O-P distance of 1.2 nm).

The electrostatic potential surfaces of HuAChE, 2-PAM, DZP and DZP_anion were calculated with the APBS software,^18,19 which may be added as a plugin to PYMOL.²⁰ The method of Poisson-Boltzmann,^19,21-23 with a dielectric constant equal to 20, was used for calculation of the electrostatic potential surfaces of 2-PAM, DZP, DZP_anion and the protein active site. The temperature was set to 310 K for all the calculations and the variation of electrostatic potential was between –1 and + 1 kT e^-1 for the active site of HuAChE and for the compounds. Finally, the softwares MOLMOL,²⁴ PYMOL²⁰ and GRACE,²⁵ were used to visualize and analyze the results of the MD simulations.

The hardware resources used in this work were a PC AMD Atlhon 2,000 MHz and a cluster of computers composed by five 2.4 GHz Pentium IV CPUs and the real time of each MD simulation was of about 10 days.

Results and Discussion

Homology modeling and active site determination

After incorporating the missing residues from the initial crystallographic structure of HuAChE,⁶ the SWISS-MODEL server^8,9 produced our complete HuAChE model where the first loop was completed by residues Pro259, Gly260, Gly261, Thr262, Gly263 and Gly264 and the second loop was completed by residues Arg493 and Asp494. The complete enzyme model is shown in Figure 3.

The Ramachandram plot for the model (see Supplementary Information) obtained from the PROCHECK site,²⁶ presented 98.7 % of the residues at the most favorable regions and only 1.3 % of the residues in the other regions. Regarding the main chain properties of the modeled loops, no considerable bad contacts, nor C_a tetrahedron distortion, nor hydrogen-bond energy problems, were found. Moreover, the average G factor, the measure of the normality degree of the protein properties, was 0.1, which is within the permitted values for homology models. Also, there were not found any distortions of the side chain torsion angles. Accordingly, and since our complete model was based on the crystallographic structure, the model was considered appropriate to study the dynamics of the compounds inside the active site of the enzyme.

The identification of the residues belonging to the active site of HuAChE was obtained by the sequence alignment between HuAChE and the crystallographic structure of AChE from Torpedo californica (TcAChE, PDB code 1EA5)²⁷ and also by comparison with the residues reported by Shafferman et al.,²⁸ Enyedy et al.²⁹ and Ariel et al.³⁰ For example, 84% of the aminoacids of the active site of HuAChE were identical to the residues of the TcAChE active site as shown in Table 2.

Molecular dynamics simulations

In order to have an insight on the behavior of the three compounds inside and outside the inhibited active site of HuAChE, there were carried out manual dockings of 2-PAM, DZP and DZP_anion inside the active site of HuAChE bound to GA, halfway inside the well that conducts to the active site, and at the entrance of the well. It is important to notice here that in a previous work, Castro and Figueroa-Villar,⁴ performed a complete conformational and docking study of 2-PAM inside TcAChE active site inhibited by GA and established the conformation used in that work as the most favorable for interaction with SGA. Accordingly, the optimal complex of that work was used as starting position for the dynamics studies inside the active site of HuAChE in the present work. The position of the compounds during the different MD simulations was monitored using the distance between the phosphorous atom of SGA and the oxygen atom of the oximes. The initial distances for the three different MD simulations for each compound were 0.30 nm, 0.50 nm and 1.26 nm. These nine HuAChE-GA-ANT systems were submitted to 50 ps of position-restrained MD (PRMD) followed by 1000 ps of MD simulations with no restriction, as described in the methodology section. Then, there were conducted calculations of temporal RMSD on each system for 1000 frames generated at each 1 ps of MD simulation. In this case the result is a unique general value for each system monitored throughout time. Taking into account that the complexes could float inside the water box, each frame was adjusted to the former by the minimum squares method when calculating the standard deviation. In Figure 4, one can see the system equilibration after 250 ps by the stabilization of the RMSD for the systems HuAChE-GA-ANT. This behavior was common to all MD simulations, with deviations never over 0.25 nm after stabilization. We also observed that, for all simulations, the temporal RMSD of HuAChE-GA-ANT practically does not change with time, thus indicating the dynamic stability of the systems.

The spatial RMSD on each aminoacid residue was also calculated in the time range of 0 to 1000 ps, at each 20 ps, totalizing 50 frames. Figure 5 shows the qualitative illustration, in the sausage representation, of the spatial RMSD for the HuAChE-GA-ANT system. Analyzing Figure 5 we can observe that the most mobile regions along the MD simulations (major values of RMSD and major thickness in the sausage representation) correspond to the residues near the N and C-teminal of each monomer and to the loops. On the other hand, the residues in the active site region and at the a-helixes and b-sheets present lower values, revealing to be more stable regions, as expected.

Analysis of the dynamic behavior of 2-PAM, DZP and DZP_anion inside and outside the active site of HuAChE

In order to analyze and compare the dynamic behavior of 2-PAM and the neutral and anionic forms of DZP, inside and outside the active site of HuAChE along the MD simulations, we extracted one frame at each 20 ps of MD simulation totalizing 50 frames for each system. The interactions of each compound with HuAChE active site were analyzed by evaluating the distance variation between the phosphorous atom of SGA and the oxime oxygen atoms of 2-PAM and the neutral and anionic forms of DZP. As described in the experimental section, this procedure afforded, for each compound, three HuAChE-GA-ANT systems further submitted to 1000 ps of MD simulations.

Our results showed that, when placed inside the active site (O-P distance 0.3 nm) 2-PAM and DZP move away from SGA, stabilizing at 0.5 nm, while DZP_anion moves further away, stabilizing at 0.7 nm (Figure 6a). When the initial position is 0.5 nm away from SGA, 2-PAM remains at that distance while DZP shows a tendency to move further apart, stabilizing at 0.7 nm and DZP_anion clearly moves away to 0.8 nm. (Figure 6b). The plot in Figure 6c show the distance variations when the compounds are placed outside the active site at 1.26 nm from SGA. It is clear that in this case 2-PAM approaches the active site moving about 0.50 nm inside the well. On the other hand, DZP and DZP_anion move away from SGA by 0.30 and 0.60 nm, respectively, showing that these molecules are being expelled from the entrance of the HuAChE active site conducing well.

When the compounds are placed the closest to SGA, all move away. The fact that even 2-PAM moves away from SGA in this case simply indicates that 0.30 nm is too close to the phosphorylated serine and that some steric hindrance repulsions are active. However, 2-PAM is able to periodically approach SGA during the MD simulation, thus showing that it could be able to react with SGA to free the serine residue. When all the compounds are at an initial 0.50 nm from SGA, 2-PAM remains close to that position while DZP and its anionic form moves away. In this case, the aromatic rings of 2-PAM, DZP and anionic DZP remain, along the MD simulations, buried in a negative pocket defined by the side chains of Asp74, Trp86, Tyr124, Ser125, Tyr337, Tyr341 and Tyr449. 2-PAM remained all the time inside the active site, as it was located at the proper position for interaction with Trp86 (Figure 7a). On the other hand, DZP was located in a favorable position for hydrogen bond formation with Tyr337, thus also remaining most of the time inside the active site (Figure 7b). Despite its great oscillations, we believe that the anionic form of DZP was not further moved away from the active site simply because it was able to establish electrostatic interactions with Tyr341 and Phe297, and a hydrogen bond with Tyr124 (Figure 7c), but the clear tendency of this compound is to eventually leave the well of access to the active site of HuAChE. In fact, when all the compounds are set just outside the well, only 2-PAM moves inside the well while DZP and DZP_anion are completely expelled.

To explain the differences in the dynamic behavior observed for 2-PAM and its deaza analogues, we decided to compute the electrostatic contours for HuAChE, 2-PAM, DZP and DZP_anion using the APBS software.^18,19 In those results, presented in Figure 8, the blue surface connects all points having a positive potential and the red surface connects all points having a negative potential while the light gray surfaces are neutral. The analysis of Figure 8 reveals that the active site cavity of HuAChE is a continuous region of negative electrostatic potential while 2-PAM presents a totally positive electrostatic potential surface, DZP_anion is almost totally negative and DZP is practically neutral. This suggests that molecules with a net positive charge like 2-PAM could be attracted to the active site cavity of HuAChE while neutral or negative molecules, like its DZP analogues, could be repelled, in perfect agreement with the results of the MD simulations. However, since some neutral molecules are able to enter the active site of this enzyme it is still possible that DZP could, by other means than simple electrostatic interactions, reach the active site of HuAChE.

The calculation of the relative binding affinities for the complexes would provide a more quantitative measure of the electrostatics, desolvation and apolar contributions of the interaction of the antidotes with the inhibited enzyme, however, in this work we used the electrostatic contours in a qualitative way, and the binding affinities were not calculated but will be discussed in detail in another article.³¹

Conclusions

The molecular dynamics simulations of 2-PAM, DZP and DZP_anion at different positions in relation to the phosphorylated Ser203 of HuAChE inhibited with the neurotoxic warfare agent tabun, showed that the positive charge of 2-PAM is very important for the transportation of this compound to the active site of this enzyme. It was shown that the neutral analogue DZP and its anion have a tendency to leave the active site of HuAChE, a process that, in this study, was interrupted by the establishment of intermolecular interactions of those compounds with some side chain residues of aminoacids in the gorge. Tai et al.³² report an opening and closing movement of the gorge in studies of MD simulations with mouse AChE. Bui et al.³³ have shown that this local conformational fluctuations of the gorge residues and large scale collective motions of the protein correlate highly with the crossing of tetramethylammonium by the AChE gorge. Those fluctuations could also be involved in the retention of DZP inside the gorge.

Despite the possibility of neutral compounds being able to reach the active site of AChE, our results show that the neutral and the anionic analogues of 2-PAM would not be able, or would be seriously hindered in relation to 2-PAM, to enter the well that conducts to the active site of HuAChE, due to their electrostatic repulsion or the lack of electrostatic attraction with the well. We believe that those results could be extrapolated to all other cationic oximes used as antidotes for organophosphorous compounds intoxication.

According to our results, when 2-PAM is initially set at 0.3 nm from the phosphorylated serine residue it moves 0.2 nm away, to what it looks to be its most stable position inside the active site of tabun-inhibited HuAChE. At this new position, the oxygen atom of 2-PAM is 0.5 nm away from the phosphorous atom of the phosphorylated serine, thus making necessary that the antidote moves to a distance of about 0.3 nm in order to start the dephosphorylation reaction. This proximity between 2-PAM and the phosphorylated serine is periodically reached after the dynamic stabilization of the system.

The results obtained in this work suggest that oximes that could change their charge from positive to neutral or negative after phosphorylation would be able to more effectively leave the active site of AChE. This behavior would minimize the reversion of the dephosphorylation of AChE, a phenomenon that is a real problem with the actual antidotes. Those results are now being used to design new and potentially more effective antidotes for organophosphate intoxication.

Acknowledgments

The authors wish to thank the Brazilian agencies CNPq, FAPERJ and CAPES (Pró-Defesa) for financial support.

Supplementary Information

The Ramachandram plot for the model obtained from the PROCHECK site is avaliable free of charge at http://jbcs.sbq.org.br, as PDF file.

12. http://biotech.ebi.ac.uk:8400/, accessed in August 2004.

31. We thank the referee for this suggestion.

Received: September 12, 2005

Published on the web: July 6, 2006

FAPESP helped in meeting the publication costs of this article.

Supplementary Information

• 1. Tafuri, J.; Roberts, J.; Ann. Emerg. Med. 1987, 16, 193.
• 2. Koelle, G.B.; Fundam. Appl. Toxicol. 1981, 1, 129.
• 3. Eddleston, M.; Szinicz, L.; Eyer P.; Buckley, N.; QJM-Ass. Int. J. Med. 2002, 95, 275.
• 4. Castro, A.T.; Figueroa-Villar, J.D.; Int. J. Quant. Chem. 2002, 89,143.
• 5. da Silva, G.R.; Júnior, I.B.; Figueroa-Villar, J.D.; Int. J. Quant. Chem. 2005, 105, 260.
• 6. Kryger, G.; Harel, M.; Giles, K.; Toker, L.; Velan, B.; Lazar, A.; Kronman, C.; Barak, D.; Ariel, N.; Shafferman, A.; Silman, I.; Sussman, J. L.; Acta Crystallogr. 2000, 56, 1385.
• 7. Guex, N.; Peitsch, M. C.; Electrophoresis 1997, 18, 2714.
• 8. Schwede, T.; Kopp, J.; Guex, N.; Peitsch, M. C.; Nucl. Ac. Res. 2003, 31, 3381.
• 9. Peitsch, M. C.; Biotechnology 1995, 13, 656.
• 10. Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R.; Comp. Phys. Comm. 1995, 91, 43.
• 11. Lindahl, E.; Hess, B.; van der Spoel, D. J.; Mol. Mod. 2001, 7, 306.
• 13. Schuettelkopf, A.W.; van Aalten, D. M. F.; Acta Crystallogr. 2004, D60, 1355.
• 14. van der Spoel, D.; van Buuren, A. R.; Apol, E.; Meulenhoff, P. J.; Tieleman, D. P.; Sijbers, A. L. T. M.; Hess, B.; Feenstra, K. A.; Lindahl, E.; van Drunen, R.; Berendsen, H. J. C.; Gromacs User Manual, version 3.1.1, Groningen, The Netherlands, 2002.
• 15. Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S. J.; Windus, T. L.; Dupuis, M.; Montgomery, J. A.; J. Comput. Chem. 1993, 14, 1347.
• 16. Darden, T.; York, D.; Pedersen, L.; J. Chem. Phys. 1993, 98, 10089.
• 17. Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G.; J. Chem. Phys. 1995, 103, 8577.
• 18. Baker, N. A.; APBS 0.3.2 User Guide Adaptive Poisson-Boltzmann Solver, Washington University in St. Louis Department of Biochemistry and Molecular Biophysics Center for Computational Biology, Washington University, St. Louis, 2004.
• 19. Baker, N. A.; Sept, D.; Joseph, S.; Holst, M. J.; McCammon, J. A.; Proc. Natl. Acad. Sci. USA 2001, 98, 10037.
• 20. DeLano, W. L.; Bromberg, S.; PyMOL User's Guide, DeLano Scientific LLC: San Francisco-CA, 2004.
• 21. Holst, M.; Baker, N.; Wang, F.; J. Comput. Chem. 2000, 21, 1319.
• 22. Baker, N.; Holst, M.; Wang, F.; J. Comput. Chem., 2000, 21, 1343.
• 23. Baker, N. A.; Sept, D.; Holst, M. J.; McCammon, J. A.; IBM J. Res. Dev. 2001, 45, 427.
• 24. Koradi, R.; Billeter, M.; Wüthrich, K.; J. Mol. Graph. 1996, 14, 51.
• 25. Turner, P. J.; Grace Development Team [http://plasma-gate.weizmann.ac.il/Grace/], Copyright (c) 1996-2004 accessed in august 2004.
• 26. Laskowski, R. A.; Macarthur, M. W.; Moss, D. S.; Thornton, J. M.; J. Appl. Cryst. 1993, 26, 283.
• 27. Koellner, G.; Kryger, G.; Millard, C. B.; Silman, I.; Sussman, J. L.; Steiner, T.; J. Mol. Biol. 2000, 296, 713.
• 28. Shafferman, A.; Ordentlich, A.; Barak, D.; Stein, D.; Ariel, N.; Velan, B.; Biochem. J. 1996, 318, 833.
• 29. Enyedy, I. J.; Kovach, I. M.; Bencsura, A.; Biochem. J. 2001, 353, 645.
• 30. Ariel, N. Ordentlich, A.; Barak, D.; Bino, T.; Velan, B.; Shafferman, A.; Biochem. J. 1998, 335, 95.
• 32. Bui, J. M.; Henchman, R. H.; McCammon, J. A.; Biophys. J. 2003, 85, 2267.
• 33. Tai, K.; Shen, T.; Börjesson, U.; Philipoulos, M.; McCammon, J. A.; Biophys. J. 2001, 81, 715.

Publication Dates

• Publication in this collection
  01 Dec 2006
• Date of issue
  Oct 2006

History

• Received
  12 Sept 2005
• Accepted
  06 July 2006