Molecular Dynamics of the Interaction of Pralidoxime and Deazapralidoxime with Acetylcholinesterase Inhibited by the Neurotoxic Agent Tabun

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,2However, 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,2e 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,2he 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. 3 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, 3 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). 3n a previous work 4 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-1methylpyridine)-1-(4-carboxyiminopyridine)-dimethyl ether hydrochloride). 5n 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., 6 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. 7Further, 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, 10 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 12 for validation.
2-PAM, DZP and DZP anion were manually fitted at a conformation proposed by Castro et al. 4 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 13 in order to obtain the topologies needed for posterior MD simulations with the GROMACS 3.2.1 package. 10,11As the GROMOS 87 force field 10 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.rtpfile) 14 in order to create a residue having a covalent bond between GA and Ser203 (aminoacids.datfile). 14Charge calculations for 2-PAM, DZP, DZP anion and SGA atoms were carried out using the software GAMESS US 15 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 3 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. 20he method of Poisson-Boltzmann, 19,[21][22][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, 24 PYMOL 20 and GRACE, 25 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.

Homology modeling and active site determination
After incorporating the missing residues from the initial crystallographic structure of HuAChE, 6 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, 26 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 α 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) 27 and also by comparison with the residues reported by Shafferman et al., 28 Enyedy et al. 29 and Ariel et al. 30 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, 4 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.   he 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 α-helixes and β-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    20 and GRACE. 25oo 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,19In 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. 31

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. 32 report an opening and closing movement of the gorge in studies of MD simulations with mouse AChE.Bui et al. 33 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.

Figure 1 .
Figure 1.Inhibition, desinhibition and aging of acetylcholinesterase. X is the leaving group.

Figure 3 .
Figure 3. Human acetylcholinesterase with the two modeled loops.Figure prepared with the program PYMOL.20

Figure 5 .
Figure 5. Sausage representation of the spatial RMSD for the system HuAChE-GA-ANT.Figure prepared with the program MOLMOL.24

Figure 6 .
Figure 6.Distance variations between the phosphorous atom of GA and the oxygen of 2-PAM (black), DZP (dark gray) and DZP anion (gray) when docked: a) Inside the HuAChE active side; b) Halfway inside the well that conducts to the active site and c) At the entrance of the well.Figures prepared with softwaresPYMOL20 and GRACE.25 Figure 6.Distance variations between the phosphorous atom of GA and the oxygen of 2-PAM (black), DZP (dark gray) and DZP anion (gray) when docked: a) Inside the HuAChE active side; b) Halfway inside the well that conducts to the active site and c) At the entrance of the well.Figures prepared with softwaresPYMOL20 and GRACE.25

Figure 7 .
Figure 7. Interactions of: a) 2-PAM, b) DZP and c) DZP anion with residues of the active site of HuAChE when docked at 0.5 nm from SGA. Figure prepared with the software PYMOL.20

Figure S2 .
Figure S2.Distance variations between the phosphorous atom of GA and the oxygen of 2-PAM, DZP and DZPanion, along the MD simulations, when docked inside the HuAChE active side; Figures prepared with the software PyMOL.

Figure S3 .
Figure S3.Distance variations between the phosphorous atom of GA and the oxygen of 2-PAM, DZP and DZPanion, along the MD simulations, when docked halfway inside the well that conducts to the active site of HuAChE; Figure prepared with the software PyMOL.

Table 2 .
Comparison between the active site's aminoacids of HuAChE and TcAChE