Docking and molecular dynamics studies of peripheral site ligand–oximes as reactivators of sarin-inhibited human acetylcholinesterase

In the present work, we performed docking and molecular dynamics simulations studies on two groups of long-tailored oximes designed as peripheral site binders of acetylcholinesterase (AChE) and potential penetrators on the blood brain barrier. Our studies permitted to determine how the tails anchor in the peripheral site of sarin-inhibited human AChE, and which aminoacids are important to their stabilization. Also the energy values obtained in the docking studies corroborated quite well with the experimental results obtained before for these oximes.


Introduction
Organophosphorus (OP) compounds are mainly used today as pesticides in agriculture and insecticides. However, some OP compounds are too toxic to be used for these purposes. They are called nerve agents, a class of warfare weapons, able to phosphorylate a serine residue in the active site of the enzyme acetylcholinesterase (AChE), blocking the hydrolysis of the neurotransmitter acetylcholine (ACh). This process leads to what is called cholinergic syndrome, whose main consequences are involuntary muscular movements and death by cardiorespiratory stop. This reaction can be reverted by oximes, which hydroxyl group carries out a nucleophilic attack on the phosphorylated serine residue, removing the phosphate moiety and reactivating AChE, as shown in Figure 1.
Despite the existence of a number of oximes able to reactivate AChE in the blood stream, the reactivation of inhibited AChE in the central nervous system (CNS) is still a challenge. This happens because the more effective oximes usually are positively charged. These oximes have a good affinity for the catalytic site of AChE (A-site), promoting a faster reactivation, but do not penetrate in the lipophilic blood-brain barrier (BBB) to protect the CNS. On the other hand, non-charged oximes have the opposite behavior. Based on this fact, Koning and co-workers have designed, synthesized, and tested 14 new oximes, being 7 charged (de Koning, Joosen, Noort, van Zuylen, &Tromp, 2011) and7 uncharged (de Koning, van Grol, &, with lipophilic tails aimed to enable the penetration in the BBB. These tails are also intended to bind to the peripheral site of AChE (P-site), known as being lipophilic. Koning and co-workers have also shown that when the affinity between the oxime and the P-site increases, the affinity for the whole system also increases, leading to stabilization. Consequently, the reactivation rate also grows up .
All the oximes studied by co-workers (2011, 2011) have shown significant reactivation rates face to sarin-inhibited human AChE (HssAChE). This result was attributed by the authors to the interactions of their tails with the P-site. However, the residues of the P-site that effectively interact with these oximes have not been described yet. The purpose of the present study is, therefore, to perform docking and molecular dynamic simulations on a model of sarin-inhibited HssAChE, using the oximes proposed by co-workers (2011, 2011), in order to determine the P-site residues interacting with their tails in order to help in the improvement of the structures and the design of more effective oximes as BBB penetrators.

Oximes
The structures of the oximes proposed by co-workers (2011, 2011) and studied in the present work, are presented in Figure 2. All of them but 4-PAM and oxime 2d were planned to have a spacer and a peripheral site ligand moiety. Oxime 4-PAM was used in the experimental studies  to compare with 2-PAM (pralidoxime), which is the most used. Oxime 2d was used to compare with the quaternary oximes 2a, 2b, and 2c .

Docking energies calculations
The model of HssAChE inhibited by sarin used in this work was constructed through mutations in the crystal structure of Mus musculus AChE (MmAChE) complexed with the oxime HI-6 available in the RCSB PDB Protein Data Bank (http://www.rcsb.org/pdb/home/home.do), under the code 2WHP, using the software SPDBViewer (Guex & Peitsch, 1997). Target and template presented 88.5% of homological identity and the residues of active sites are 100% conserved. The 3D structures of each compound in Figure 2 were built using the program PC Spartan 08 ® (Hehre, Deppmeier, & Klunzinger, 1999) and their partial atomic charges calculated by the RM1 (Recife Model 1) semi-empirical method (Rocha, Freire, Simas, & Stewart, 2006). After charges calculations, the docking energies of the ligands bound to the active site of HssAChE were obtained using the software Molegro Virtual Docker (MVD) ® (Thomsen & Christensen, 2006) using the algoritm MolDock Score, an adaptation of the algorithm Differential Evolution (Thomsen & Christensen, 2006). The MolDock Score energy, E score , is defined by Equation (1), where E inter is the ligand-protein interaction energy, and E intra is the internal energy of the ligand. E inter is calculated according to Equation (2).
The E PLP term is a "piecewise linear potential" which uses two different parameters, one for the approximation of the steric term (van der Waals) between atoms, and another, which is a potential for hydrogen bonds. It describes the electrostatic interactions between charged atoms (Thomsen & Christensen, 2006). E intra is calculated according to Equation (3).
The first term in Equation (3) calculates all the energies involving pairs of atoms of the ligand, except those connected by two bonds. The second term refers to the torsional energy, where θ is the torsional angle of the bond. The average of the torsional energy bond contributions is used if several torsions could be determined. The last term, E clash , assigns a penalty of 1000 kcal mol −1 if the distance between two heavy atoms (more than two bonds apart) is smaller than 2.0 Å, ignoring infeasible ligand conformations (Thomsen & Christensen, 2006). The energies of Equations (1)-(3) are given in kcal mol −1 . The docking accuracy of MolDock has been evaluated by docking flexible ligands to 77 protein targets. Mol-Dock was able to identify the correct binding mode of 87% of the complexes. In comparison, the accuracy of docking programs like Glide and Surflex is 82 and 75%, respectively. FlexX obtained 58% and GOLD 78% on subsets containing 76 and 55 cases, respectively (Thomsen & Christensen, 2006). The compounds were docked in the HssAChE binding site using MVD ® but first, a re-docking procedure was performed using the structure of oxime HI-6 from the crystal in order to validate the protocol used. The binding site was restricted into a sphere with a radius of 15 Å and the residues within a radius of 10 Å were considered flexible. Due to the stochastic nature of the docking algorithm, about 15 runs were performed for each compound with 30 poses (conformation and orientation of the ligand) returned to the analysis of the distance between the P atom of sarin and the O atom of the oxime. The best conformation of each compound was selected according to this distance, the energy and the interactions observed. These conformations were used in further steps of molecular dynamics (MD) simulations.

MD simulation
Before performing the MD simulations it was necessary to parameterize the ligands, so that they could be recognized by the force field OPLS/AA (Jorgensen, Maxwell, & Tirado-Rives, 1996) from the program GROMACS 4.0 (Hess, Kutzner, Van Der Spoel, & Lindahl, 2008). To obtain the parameters and topologies for the referred compounds, it was used AnteChamber PYthon Parcer InterfacE (ACPYPE) (Da Silva & Vranken, 2012) in association with MKTOP (Ribeiro, Horta, & de Alencastro, 2008). ACPYPE is a tool based on PyThon programming language to use ANTECHAM-BER (currently bundled in AmberTools version 1.4) (Wang, Wang, Kollman, & Case, 2006) to generate parameters and topologies for chemical compounds and to interface with other Python applications like CCPN tools (Vranken et al., 2005) or ARIA2 (Rieping et al., 2007). ACPYPE is currently able to generate output files for the following MM softwares: CNS/XPLOR (Brunger, 2007;Schwieters, Kuszewski, & Clore, 2006), GRO-MACS (Hess et al., 2008), CHARMM (Brooks et al., 1983), and AMBER (Hornak et al., 2006). OPLS/AA (Jorgensen et al., 1996) force field parameters for the ligands were used and the atomic partial charges were calculated by semi-empirical quantum chemistry SQM (Walker, Crowley, & Case, 2008) (by ANTECHAM-BER) according to AM1-BCC parameters. SQM was modified to include 6 decimals of precision instead of the default 3. The enzyme/Sarin/ligands complexes were simulated using the GROMACS 4.0 (Hess et al., 2008) package, in a cubic box (924.67 nm 3 ) containing approximated 28,000 TiP3P water molecules with periodic boundary conditions, and this was the starting configuration. The minimization algorithms used were steepest descent with position restrained (PR) of ligands and protein, in order to accommodate water molecules to avoid overlapping of van der Waals Radis, with a convergence criterion of 100.00 Kcal mol −1 Å −1 , steepest descent without PR to flexibly the system, and conjugate gradients to look for the local minimum. The minimized complexes were then submitted to MD simulations in two steps. Initially, we performed 500 ps of MD at 310 K with PR for the entire system, except the water molecules, in order to ensure a balance of the solvent molecules around the residues of the enzyme. Subsequently, 10,000 ps of MD were performed at 310 K without any restriction, using 2 fs of integration time and a cut-off of 10 Å for long-distance interactions. Counter ions were added to the systems in order to neutralize the total charge of the enzyme-ligand complexes.
In order to analyze the structures generated after the optimization and MD steps, we used the VMD program (Humphrey, Dalke, & Schulten, 1996). Plots of variation of total energy, distance, the variation of random mean square deviation (RMSD), and H-Bonds formed along the MD simulation were generated with the Origin program (Edwards, 2002). Qualitative spatial RMSD pictures were generated in the MolMol (Koradi, Billeter, & Wüthrich, 1996) and the figures of the frames of MD simulations were generated in the PyMOL program (Warren, 2002).

Free energies calculations
In order to check our conclusions about the binding energy of the compounds docked with the enzymes, we used the tool g_mmpbsa which use the molecular mechanics Poisson−Boltzmann surface area method (MM-PBSA) (Kumari, Kumar, & Lynn, 2014). The auto correlation time of the system HssAChE/sarin/oximes was about 450 ps (calculate with g_analyze), so we sampled the trajectory each 500 ps where each frame is not correlated with the previous.
MM-PBSA as implemented by Kumari et al. (2014) is designed to analyze the solvation properties of small and macro-molecules such as proteins, nucleic acids, and other complex systems. This method combine three contributions to the free energy: the first one correspond to the potential energy in vacuum and includes bonded terms (such as bond, angle, and torsion energies) as well as non-bonded terms (such as van der Waals and electrostatic interactions). The solvation effects are considered in the second term, which includes the sum of two energy terms, i.e. the polar and non-polar solvation energies, using an implicit solvation model. The last term includes the entropic effect with the complex in the gas phase (Kumari et al., 2014).
In order to quantify and determine the relevance of each amino acid in the interactions protein-ligand, we used the intermolecular surface contact (ISC) (Cuya Guizado, Louro, & Anteneodo, 2012). The ISC could furnish information on the interactions between the oximes and given HssAche residues. The ISC was obtained by means of the SURF software (Connolly, 1983). This software is based on the Connolly algorithm, using a 1.4 Å test radius and 3 point/Å. The ISC is determined by the intersection of the surfaces accessible to the solvent of both the protein and the ligand.
3. Results and discussion 3.1. Docking studies Figure 3 shows the re-docking result for the best docked structure of HI-6 inside HssAChE. RMSD for the superposition of the non-hydrogen atoms was 1.230 Å. Keeping in mind that a RMSD value under 2.0 Å is considered acceptable (Kontoyianni, McClellan, & Sokol, 2004), this result validates the docking protocol used.
The ligand-protein interaction energies (intermolecular and H-bonds) were calculated in order to get a better understanding of the variations between the binding modes of each compound and the molecular factors responsible for the reactivation. Table 1 lists the residues present in these interactions, the energy values, the possible interaction residues in the P-site area, interatomic P-O distances, and the experimental results of reactivation for HssAChE obtained by . Analyzing the plots in Figure 4, it can be noticed a very good correlation between the theoretical results obtained and experimental data.
The ligands with the smallest intermolecular energy values (oxime 1c, for charged ones, and oxime 2b, for non-charged ones), considered experimentally as the most promising for HssAChE, correspond to the highest reactivation percentages. It points to great stabilities in active site. Table 1 also shows that the P-O distance tends to reduce as long as energy is smaller, with few exceptions. This proximity is very significant to the reactivation process and also points to corroboration with experimental data. As it was concluded by , the hydrophobic interactions in P-site area are very important to increase the affinity of the ligand with the enzyme and, consequently, the reactivation potential. The docking studies appear to confirm that. The oximes with aromatic rings in P-site area (1a-c and 2a-c), that presented better experimental results, also presented smaller intermolecular energy values and smaller P-O distance in general than chlorine oximes (3a-c and 5a-c). Figure 5 shows the best docking conformations for oximes 1c and 3c (charged), respectively. Figure 6 shows the best docking conformation for oximes 2b and 5b (non-charged). Some estimated P-site interaction residues are the same in almost all the dockings experiments, especially between the ligands that are different from the other just in the length of the spacer. And it is logical because they have the same aromatic rings in the P-site area. So, it can be inferred that Trp86, Tyr124, Trp286, Phe297, Tyr337, and Tyr341 are the most contributive residues in hydrophobic interaction and ligands stabilization in the P-site area.

MD simulation
After docking, the oximes were submitted to 10 ns of MD simulations in order to check their stability inside HssAChE and verify which residues interact with the Psite of HssAChE.
Results of the MD simulations showed that for all systems, the total energy tends to be stable after 1.0 ns of simulation. This behavior points to constant energy and suggests structural stabilization as shown in Supplementary Figure S1.
The temporal RMSD was calculated on all atoms of each complex to 500 frames at every 20 ps, during the 10 ns of simulation. Considering that the complexes could fluctuate in the box, each frame was adjusted by the least squares method to its previous one for the calculation of the standard deviation. In Figure 7, it is possible to observe the equilibration of the simulation for the systems HssAChE/oximes (1a-c) and HssAChE/oximes (2a-c) around the initial 2000 ps. This behavior was common to all simulations, with deviations never     exceeding .22 nm (2.2 Å) and .10 nm (1.0 Å) for protein and ligand, respectively. This result suggests that the compounds accommodate well inside the active site along the 10 ns of simulation, showing stabilization of the system and confirming the results obtained by the total energy calculations previously described.
The root mean square fluctuations (RMSF) on each amino acid were also calculated and are illustrated in Figure 8, where the red boxes contain the amino acids in the catalytic and P-sites. For these sites all the RMSF profiles were similar except for oxime 2b, where deviation was of about 2 Å, shown by the difference in size of the peak of this oxime related to the others around residue 250, in Figure 8. The residues responsible for this increasing of flexibility are Gly252, Cys253, Pro254, Pro255, Gly256, Thr258, Glu259, and Gly260.
In order to facilitate the visualization of the dynamic behavior of each oxime within the active site and P-site, we evaluated the variation of the distances of the residues directly involved in interactions with each oxime for each system, according to the results obtained in the docking studies. In general, the main interactions observed in the docking studies were also seen in the MD simulations. Considering this, we will discuss in more detail only two charged oximes (one with aromatic rings in P-site area and another without it) and two noncharged oximes (one with aromatic rings in P-site area and another without it).

Oxime 1a
About three hydrogen bonds were formed with oxime 1a during the 10 ns of MD simulation (Supplementary Figure S2) and Glu202, Tyr337, and Tyr341 were the main interacting residues. The average distance between the mass center of the oxime and these residues were 1.23, .71, and .62 nm, respectively. The distance between the mass center of the oxime and the P-site interaction residues pointed in docking studies (Trp86 and Trp286) were .92 and .99 nm (Supplementary Figure S3). The superposition of frames for oxime 1a and the interacting residues after 10 ns of MD simulation is shown in Figure 9.

Oxime 3a
Only one hydrogen bond in average was formed with oxime 3a during the 10 ns of MD simulation (Supplementary Figure S2). Tyr124 and Ser298 were pointed in docking studies as interaction residues. The hydrogen bond with Tyr124 was observed during the MD simulation, but it did not last up to the end of the simulation.   With Ser298, no hydrogen bond was observed. This suggests a low stability of this oxime inside HssAChE. The average distance between the mass center of the oxime and the main interaction residues of active site and P-site showed a low stability. Except for Phe297 all the residues were far from the oxime. The variation of the interatomic distance between the sarin's phosphorus and the oxime's oxygen showed that these two atoms were kept far from each other during all the simulation time (Supplementary Figure S4). These results were already expected because oxime 3a was pointed in experimental studies as a poor reactivator. The superposition of frames for oxime 3a and the interacting residues after 10 ns of MD simulation is shown in Figure 10.

Oxime 2a
Oxime 2a formed up to 3 hydrogen bonds during the 10 ns of MD simulation (Supplementary Figure S2), but only the bond with Glu202 was maintained. The average distances between the mass center of this oxime and residues Tyr124, Phe295, and Phe297 (Supplementary Figure S5) were kept stable and close to the oxime during the simulation time, pointing to the stability of this oxime in the P-site. Besides the sarin's phosphorus and oxime's oxygen interatomic distance was kept stable and close to 5 Å during the simulation time. These suggest that oxime 2a is a promising reactivator, as shown in experimental studies. The superposition of frames for oxime 2a and the interacting residues after 10 ns of MD simulation is shown in Figure 11.

Oxime 5a
Oxime 5a formed about four hydrogen bonds during 10 ns of MD simulation (Supplementary Figure S2), but only the bond with Gly120 remained till the end of the simulation time. The average distances from the mass center of oxime 5a to Trp86, Gly120, and Ser125 were close, but not stable (Supplementary Figure S6). Distances to Trp286 and Tyr337 showed some stability, but they were far from the oxime. Sarin's phosphorus and oxime's oxygen interatomic distance increased during the simulation (Supplementary Figure S4). All these facts suggest that oxime 5a is a poor reactivator, corroborating the experimental results. The superposition of frames for Figure 11. Superposition of frames of oxime 2a and the interacting residues after 10 ns of MD simulation.  oxime 5a and the interacting residues after 10 ns of MD simulation is shown in Figure 12.

Oximes binding energies
In order to validate the docking energies we used the molecular mechanics Poisson-Boltzman method surface area (MM-PBSA) (Kumari et al., 2014) to calculate the binding energies of the oximes inside HssAChE/sarin. Our results showed correlation with the docking energies and the experimental % of reactivation of the oximes.
As can be seen in Figures 13 and 14 the linear fit showed correlations with R 2 = .89 and .72, respectively.
Oximes 1c, 2a, and 2b are the best candidates in terms of energy and % of reactivation. In the case of oximes 2a and 2b the differences in the docking energy is about .2 kcal mol −1 . It is insufficient to evaluate the best one using only molecular docking. However, the MM-PBSA method can distinguish and select the best candidate. As expected the MM-PBSA method confirmed our docking protocol and methodology and pointed to the best available candidates.

Intermolecular surface contact
The time-averaged ISC was computed for the trajectories within the 5-10 ns interval of MD simulation. The   amino acids with ISC > 20 Å 2 for the most promissory oximes 1c and 2b are displayed in Figure 15. All of them belong to the catalytic site. ISC was also calculated for all the oximes as shown in Figure 16. Our analysis suggests that the amino acids interacting more strongly with the oximes, in decreasing order of ISC, are: Trp82, Tyr120, Trp282, Glu288, Tyr337, His283, and Phe334.
Docking and MD simulations showed that the tryptophanes and tyrosines play an important role in the oxime binding. Additionally, MD based on rigorous statistical sampling revealed, through ISC, the residues in close contact with the oximes. Figure 16 shows the correlations between ISC and % of reactivation for the oximes studied. The upper region shows that the most promissory inhibitors in terms of energy and % of reactivation are, also, the oximes with larger values of ISC while the lower region concentrates the less promissory reactivators and the lower values of ISC. This suggests that the ISC can be directly correlated to the potential of reactivation and, therefore, should be considered in the design of new oximes.

Conclusions
Our methodology seems to show a good correlationship with the experimental data reported by  and was also able to predict the residues of the P-site important for the stabilization of the oxime's tails. Oximes with a hydrophobic tail interacting with the P-site showed to be the best reactivators in docking and MD simulations, even when uncharged, fact that is good for BBB penetration. The most promising reactivators, according to our studies, were oximes 1c and 2a, corroborating with experimental data. Also, the MM-PBSA method validated our docking protocol and methodology and the ISC studies proved to be important in the design of new reactivators.

Supplementary material
The supplementary material for this paper is available online at http://dx.doi.org/10.1080/07391102.2015.1124807.