Exploration of the chlorpyrifos escape pathway from acylpeptide hydrolases using steered molecular dynamics simulations

Acylpeptide hydrolases (APH) catalyze the removal of an N-acylated amino acid from blocked peptides. APH is significantly more sensitive than acetylcholinesterase, a target of Alzheimer’s disease, to inhibition by organophosphorus (OP) compounds. Thus, OP compounds can be used as a tool to probe the physiological functions of APH. Here, we report the results of a computational study of molecular dynamics simulations of APH bound to the OP compounds and an exploration of the chlorpyrifos escape pathway using steered molecular dynamics (SMD) simulations. In addition, we apply SMD simulations to identify potential escape routes of chlorpyrifos from hydrolase hydrophobic cavities in the APH-inhibitor complex. Two previously proposed APH pathways were reliably identified by CAVER 3.0, with the estimated relative importance of P1 > P2 for its size. We identify the major pathway, P2, using SMD simulations, and Arg526, Glu88, Gly86, and Asn65 are identified as important residues for the ligand leaving via P2. These results may help in the design of APH-targeting drugs with improved efficacy, as well as in understanding APH selectivity of the inhibitor binding in the prolyl oligopeptidase family.


Introduction
Acylpeptide hydrolase (APH) belongs to the prolyl oligopeptidase (POP; EC 3.4.21.26) family of serine proteases (Barrett & Rawlings, 1992;Polgár, 2002b;Rawlings, Polgar, & Barrett, 1991). The POP family is a new group of serine peptidases and different from the classic serine proteases, trypsin (Rajapakse, Ogiwara, & Takahashi, 2014) and subtilisin (Fernández, Daleo, & Guevara, 2015), in several structural features and catalytic behaviors (Polgar, 1994(Polgar, , 2002aRennex, Hemmings, Hofsteenge, & Stone, 1991;Szeltner et al., 2003). They regulate the activity of biologically active peptides and peptide hormones, and they are implicated in diseases, including amnesia, depression, diabetes, and trypanosomiasis (Rea & Fülöp, 2006). Using the MEROPS (http://merops.sanger.ac.uk) database, which contains information on peptidase sequences, structures, substrates, and inhibitors (Rawlings, Tolle, & Barrett, 2004), the enzymes of the POP family can be classed into clan SC, family S9, and the members of which all have a catalytic triad comprising the Asp-His-Ser triad. The Asp-His-Ser catalytic triad is found in at least four structurally different classes of serine peptidase, indicating that independent evolutionary paths converged on this catalytic system (Dodson & Wlodawer, 1998). Despite these different structural frameworks, the position and geometry of the catalytic residues is similar, indicating a similar catalytic mechanism throughout the POP family.
The use of acetylcholinesterase (AChE) inhibitors in organophosphorus (OP) compounds for treating human diseases such as cognitive defects in Alzheimer's disease, has long been established (Richards, Johnson, & Ray, 2000). In principle, all serine hydrolases have the capacity to react with OP compounds, thus the characterization of members of this class of enzymes containing the POP family protein in biological systems would provide a useful resource for the identification of potential OP targets.
As mentioned above, despite the progress made in experimental and molecular simulation studies on APH and its complexes, several questions remain open. How do the OP compounds enter and leave the active site gorge of APH? What role do the residues play in the binding and unbinding of OP compounds? (Kalyaanamoorthy & Chen, 2014;Li et al., 2012;Marzinek et al., 2014;Nicolini, Frezzato, Gellini, Bizzarri, & Chelli, 2013;Vashisth & Abrams, 2008;Wang & Duan, 2009). This method requires knowledge of the unbinding pathway, which is not trivial for those proteins whose binding pocket is buried into the protein core with no obvious channels for ligand unbinding. APH is a typical example of such a protein. Due to the potential application of this method to the structural modification of chlorpyrifos (see Figure 1(a) and Figure S1), one of the OP compounds, it is important to identify the unbinding pathways through which chlorpyrifos escapes from the APH ligand-binding domain.
In this study, we will reveal the unbinding pathway and the associated unbinding mechanism of chlorpyrifos from APH. CAVER 3.0 1 (Chovancova et al., 2012) was used to identify and estimate all previously published APH tunnels for a 50-ns molecular dynamics (MD) simulation. Molecular docking was then applied to the systems of chlorpyrifos in an APH complex, and multiple SMD simulations were employed to explore the possible ligand unbinding routes. By analyzing the successful unbinding trajectories and using SMD simulations, the most likely unbinding pathway was determined, and the associated unbinding mechanism was finally elucidated.

System preparation
The initial structure of APH in the MD simulations was taken from the Protein Data Bank (PDB Id 1VE7 at 2.4 Å resolution) (Bartlam et al., 2004). The protonation states of the ionizable residues and histidines were determined from H ++ (Anandakrishnan, Aguilar, & Onufriev, 2012). Based on the calculated pKa values and the microenvironment, His90, His285, His 299, His391, His508, and His520 were set to be fully protonated at both nitrogen atoms. The initial model of the chlorpyrifos was extracted from the ChemSpider database. Geometric optimization of the chlorpyrifos was conducted at the B3LYP/6-31G** level using Gaussian 09 (Lomas, 2014;Premkumar, Jawahar, Mathavan, Kumara Dhas, & Milton Franklin Benial, 2015).

Molecular docking
CDOCKER from the Discovery Studio 2.5 software package was used for molecular docking to the APH. CDOCKER is a grid-based molecular docking method that employs CHARMM (Abuo-Rahma, Abdel-Aziz, Farag, & Kaoud, 2014). The receptor was held rigid while the substrate was allowed to flex during structural refinement. Because the crystal structure of the hydrolase had been obtained through X-ray crystallography, prior knowledge of its binding site had been acquired. Hence, it was possible to specify the placement of the substrate (inhibitor) at the active site using a binding site sphere with a radius of 8 Å. Random ligand conformations were generated from the initial ligand structure through high-temperature MD at 1000 K, followed by random rotations. The random conformations were refined by grid-based simulated annealing and a final full force field minimization. The top 10 poses were saved for comparison and analysis.

Conventional MD simulations
All of the complex systems (APH, APH-chlorpyrifos) were subjected to MD simulations with periodic boundary conditions using the GROMACS 4.5.2 software package with the SPC water model (Hess & van der Vegt, 2006). For the ligands, the Dundee PRODRG (Schüttelkopf & van Aalten, 2004) server was used to build a Gromacs topology for chlorpyrifos, and the .itp file containing the molecular information was obtained. We added the .itp file to the protein .top file, and then the protein-ligands complex simulations were again carried out using the GROMOS force field force parameter set 53A6 (Cao, Lin, Wang, & Liu, 2009). To keep each system in an electrically neutral state, the appropriate number of sodium ions was added to randomly replace the water molecules. Initially, the energies of the complex systems were relaxed through steepest descent energy minimization to exclude steric clashes or incorrect geometry. Afterwards, NVT and NPT were alternately run with position restraints on the protein and ligand to allow for relaxation of the solvent molecules in two phases. The solvent molecules were equilibrated with the fixed protein at 300 K, taking the initial velocities from a Maxwellian distribution. Subsequently, the protein and inhibitor were relaxed step by step and heated up to 358 K. The long-range electrostatics was described with the particle mesh Ewald algorithm (Darden, York, & Pedersen, 1993) with an interpolation order of 4, a grid spacing of .16 nm, and a Coulomb cut-off distance of 1.0 nm. The LINCS algorithm was used to constrain all bonds. The temperature and pressure coupling types were set to V-rescale and Parrinello-Rahman, respectively (Berendsen, Postma, van Gunsteren, DiNola, & Haak, 1984). In the NVT ensemble, the temperature of the systems reached a plateau at the desired value (reference temperature = 358 K; time constant = .1 ps). In addition, the equilibration of pressure (reference pressure = 1.0 bar; time constant = 2.0 ps) was performed using the NPT ensemble. Finally, the 50-ns MD simulations, for collecting data with a time step of 2 fs and coordinates saved every 2 ps, were initiated (RMSD plot are shown in Figure S2).
Pathways identified with CAVER 3.0 CAVER 3.0 outperforms existing software for geometry-based analysis of pathways in MD simulations (Chovancova et al., 2012). It provides detailed characteristics of individual transport pathways and their time evolution, and enables identification of pathways invisible in a static structure and investigation of the structural basis of the pathway gating mechanism. The pathways are identified with CAVER 3.0 as the paths in a graph composed of Voronoi vertices and Voronoi edges (Aurenhammer, 1991). Several pathways with nearly identical axes may be identified within one snapshot. To remove the redundant pathways, the following iterative procedure is employed. The lowest cost pathway P is selected and all pathways within a user specified distance from P are discarded. The procedure is repeated with the next remaining lowest cost pathway, until all pathways have been either selected or discarded.

SMD simulations and umbrella sampling calculations
To investigate the binding process of OP compounds to APH, an SMD simulation was carried out to pull the inhibitor, chlorpyrifos, from the active site. The SMD simulation method is a computational approach used for atomic force microscope, optical tweezer, biomembrane force probe, and surface force apparatus experiments (Mousavi, Amjad-Iranagh, Nademi, & Modarress, 2013;Zhang, Zheng, Li, & Zhang, 2013). In SMD simulations, an external force is applied to an atom (or group of atoms) along a reaction path. The SMD atom is attached to a dummy atom via a virtual spring while a reference atom (or group of atoms) is kept at a fixed location. The dummy atom can then be moved at a constant velocity along the reaction coordinate, dragging the SMD atom (Miño, Baez, & Gutierrez, 2013). The two kinds of simulations correspond to constant force SMD and constant velocity SMD .
The SMD pulling directions in the APHchlorpyrifosmethyl oxon complex were determined in the tunnel study by CAVOR 3.0 (Chovancova et al., 2012). The pulling directions were set by defining two atom groups: the Cα atom of Glu88 and its opposite residue Ala557 at the tunnel entrance (P2), and Ser157 and Leu251 surrounding the tunnel entrance (P1) (Figure 4(b)).
In the present study, SMD simulations with constant velocity were performed on the system when it reached an equilibrium state after the 50-ns MD simulation. The SMD pulling direction in the APH-chlorpyrifos complex was determined by the statistical results of the MD simulations. A spring constant of 1000 kJ mol −1 nm −2 was applied to the center of mass of the ligand for the two pathways. The P2 length is about 2 nm and the largest length of the ligand is less than .8 nm. We therefore pulled 3 nm for the side passageway, and the pulling velocity was set to .0003 nm ps −1 , which had to be slow in order to pull the ligand out successfully. For, P1 we pulled 5 nm out as the tunnel length is around 4 nm, and the pulling velocity was set to .0005 nm ps −1 . SMD simulations were carried out with the software package GROMACS 4.5.2 (Lomas, 2014), starting from the snapshot structure at 10 ns for the two passageways.
The free energy change along a reaction coordinate is termed the potential of mean force (PMF). PMF is a potential that is obtained by integrating the mean force from an ensemble of configurations. To energetically investigate the dissociation of chlorpyrifos from the active site of APH, a reaction coordinate was chosen as the distance between the mass centers of chlorpyrifos and the fixed atoms in SMD simulation. Chlorpyrifos was docked at APH active site and will be out using pulling force with proper spring constant (1000 kJ mol −1 nm −2 ) and pulling velocity. Along the reaction axis, umbrella sampling (US) windows were set at intervals of roughly .1 nm. A 3-nm reaction coordinate for P2 is separated into 30 windows for each US window. An NPT equilibration was performed followed by US production, with either procedure lasting for 500 ps, and full sampling was carried out for 10 ns for each window. As well as the P1, we separated 50 windows along the coordinate of P1. To further increase the efficiency of the calculation, PMF was calculated by averaging the outcomes of the windows independent runs.

Molecular mechanical/GBSA calculations
Ten nano second MD was used for two protein-inhibitor complex analysis. The 2000 snapshots isolated from the final 4000 ps MD trajectory with protein-substrate complex were used for the binding free energy calculation using the MM-PBSA method encoded in the Amber 10 program (Case et al., 2008). The calculation of binding free energy was carried out by means of the Molecular mechanical (MM)/GBSA method (Case et al., 2008). The method combines the MM energies with the continuum solvent approaches by the programs in the Amber 10 program, where the Sander module is used for the calculation of the MM energies, the GBSA module (Kollman, Massova, Reyes, Kuhn, & Huo, 2000) is used for the calculation of the polar solvent energy (the electrostatic contribution), and the Molsurf program (Sitkoff, Sharp, & Honig, 1994) is used for the calculation of the nonpolar solvent energy (the nonpolar contribution).

Results and discussion
The binding pose of chlorpyrifos of APH In our research, molecular docking was used for generating the enzyme-inhibitor complex structure. The three-dimensional structure (3D) of chlorpyrifos was downloaded from the ChemSpider database and the B3LYP/6-31G** level was used with Gaussian 09 (Lomas, 2014;Premkumar et al., 2015). Figure 2(b) and (c) shows the HOMO and LUMO orbits of chlorpyrifos using the Multiwfn program (Lu & Chen, 2012). It can be concluded that the active part of chlorpyrifos is the benzene group. To facilitate analysis of the electron density properties, the contour maps from the charge density deformation program (Figure 3(a) and (c)) of chlorpyrifos and p-nitrophenyl caprylate, a substrate of APH, between the adsorption structure and the noninteraction fragments (in the same position) and the interatom surface (Figure 3(b) and (d)) were calculated from the molecular orbitals using the Multiwfn program (Kollman et al., 2000). Further analysis of the electron density and charge density deformation revealed strong chemical interaction as is visualized in Figure 3. It is generally believed that for covalent bonding due to the overlap of molecular orbitals, electrons accumulate at the center of the bond (Han et al., 2013). The electron density increases in the intermediate region indicated by black solid lines in Figure 3(a) and (c). Moreover, the charge depletion between the P of chlorpyrifos and the C of p-nitrophenyl caprylate is shown in Figure 3(a) and (c) as dashed lines. All the above, reflect the weak interaction between the P of chlorpyrifos than the C of p-nitrophenyl caprylate. Hence, the P of chlorpyrifos can be more easily attacked by the nucleophilic, Ser445.
It was reported that chlorpyrifos may affect other neurotransmitters, enzymes, and cell signaling pathways, potentially at doses below those that substantially inhibit AChE and other serine hydrolase via hydrolysis (Connors et al., 2008). Chlorpyrifos functions as a substrate and can be bind strongly to APH and impossible to remove. Hence, the chlorpyrifos located at the active site pocket contained Ser445, His556, and Asp524. The docked results from CDOCKER presented 10 top poses of the chlorpyrifos in the active site (Abuo-Rahma et al., 2014). The lowest energy between the chlorpyrifos and APH was selected for further study. As seen in Figure 3, Phe371, Gly368, Ser445, Gly369, Met477, Val471, Phe485, Pro370, Trp474, Arg526, and Ile489 were important residues for chlorpyrifos binding. Asp524 was not involved in chlorpyrifos binding. The reason for this may lie in Asp524 acting as catalytic base and therefore not being involved in substrate binding. The distance between the O atom of Ser445 and the P atom of chlorpyrifos was .285 nm, and this is useful for attacks by Ser445. It was reported that Ser445, Gly369, and Arg526 are important for substrate binding (Bartlam et al., 2004;Wang, Yang, Liu, & Feng, 2006). Our results are consistent with the experimental data (Bartlam et al., 2004;Wang et al., 2006). From Figure 6, it can be seen that Gly368, Ser445, Gly369, and Pro370 have electronic contact with chlorpyrifos, whereas Phe371, Met477, Val471, Phe485, His556, and Trp474 have strong van der Waals (VDW) contact with chlorpyrifos. In particular, there is the cation-π interaction which is recognized as an important noncovalent binding interaction relevant to structural biology (Gallivan & Dougherty, 1999), between Arg526 and chlorpyrifos. Thus, the lowest energy calculated by CDOCKER (Abuo-Rahma et al., 2014) can be assumed as the optimal binding mode for further study. The enzymeinhibitor complex was used for 50-ns MD simulation. The average 3D structure of enzyme-inhibitor complex during 50 ns was got. Seen from Figure S3 it can be seen that Arg526, His556, Phe485, Phe488, Pro370, Phe371, Ala372, Gly369, Tyr444, Gly368, and Ser445 are important residues for chlorpyrifos binding. Figure S4 indicated the distance between the O atom of Ser445 and the P atom of chlorpyrifos during 50-ns MD simulations, which floated from 0.2-0.4 nm, and this is useful for attacks by Ser445. Hence, the average 3D structure of enzyme-inhibitor complex can be used for SMD study.

Identification of pathways in MD trajectories
CAVER 3.0 was used for the analysis of 2500 snapshots from a 50-ns MD simulation of APH. In each snapshot, all possible pathways with bottleneck radius equal to or larger than .1 nm were identified, leading to a set of nearly 14,389 pathways. All nine tunnels were identified among the top ranked pathway clusters (Figure 4(a), Table 1). It was reported the main entrance of APH to the active site was via a tunnel in the propeller domain, whose diameter is approximately .7 nm by Bartlam et al. (2004). A second smaller side opening also provides access to the active site and is located between blades 1 and 2. The .6 nm wide opening is lined by the residues Asn65, Arg81, Asp82, Glu88, Asp553 Ala557, Ile558, Asn559, and Asn563 by Bartlam et al. (2004). We found a good agreement between the results of CAVER 3.0 and the experimental data (Bartlam et al., 2004): (i) the two previously proposed APH pathways were reliably identified by CAVER 3.0, with an estimated relative importance of P1 > P2 for its size; (ii) the P1 tunnel was found to be the dominant transport pathwayit was the most frequently identified collective pathway for its greater size, as it had the highest maximum and mean radii of bottlenecks by far (Table 1, Figure 4(a)); (iii) based on all of the characteristics examined, the P2 tunnel was found to be the second most important.

Analysis of bottleneck dynamics
Similarly to other tunnel characteristics, the analysis of the bottlenecks which represent promising hotspots for the modification of tunnel properties, since their substitutions have potentially the most pronounced impact on the tunnel geometry is more informative when the dynamics of the protein is considered. The list of the bottleneck residues obtained by the analysis of the MD trajectory using CAVER 3.0 reveals the structural details of the tunnel gating ( Table 2).
The most frequent bottleneck in the APH P1 tunnel is formed mainly by Ser157 (93.13% of P1 tunnel), Arg287 (92.84%), Asp158 (92.64%),and Trp250 (91.84%) residues ( Table 2). The most frequent bottleneck in the APH P2 tunnel is formed mainly by Ala557 (91.51% of P2 tunnel), Asn65 (76.62%), Asp82 (69.70%), His90 (69.23%), Glu88 (66.95%), His556 (60.98%), Gly86 (601.12%), Asn559 (57.99%), Gly555 (51.30%), Ile558 (44.76%), and Arg526 (35.66%). The proposed structural basis of gating in the bottleneck is in agreement with the results from X-ray experimental analysis, which show that the P2 tunnel is lined by the residues Asn65, Asp82, Glu88, Ala557, and Ile558 (Bartlam et al., 2004). The importance of the bottleneck residues identified by CAVER 3.0 for the catalytic properties of APH was previously demonstrated Figure 3. The important residues around chlorpyrifos calculated by Discovery studio 3.5 Client. Color purple represents for electrostatic contacts, color green represents for van der Waals contacts, and color gray represents for cation-π interactions. The top ranked collective pathways identified throughout the molecular dynamics simulation of APH by CAVER 3.0 are all depicted in one frame as pathway centerlines. The P1 and P2 tunnels were initially identified as one collective pathway using the clustering threshold of 3.5 D. A random subsample of identified tunnels P1 and P2 is shown for clarity. (b) Bottleneck dynamics and structural basis of gating in the P1 and P2 tunnel of APH. experimentally (Aurenhammer, 1991;Chovancova et al., 2012). The independently directed evolution experiment provided APH variant R526 V with 150-fold (Wang et al., 2006) higher activity towards p-nitrophenyl caprylate. With the E88A/R526K variant, both a peptidase substrate and an esterase substrate were almost abolished (Yang et al., 2009). Both of these variants carry mutations in the position of Arg526, which forms the bottleneck of the P2 tunnel. Our results are consistent with experimental data. Therefore, the study of dynamical systems enables the identification of all important bottlenecks and provides invaluable information about their relative importance.

PMF calculations
In our research, SMD simulations were performed to drive the dissociation of chlorpyrifos from the enzyme, which is the reverse process of binding. Figure 5(a) and (b) shows the changing profile of the external force vs. the simulation time for P2 and P1. As shown in the Figure 5(a), the pulling force exerted on chlorpyrifos through P2 increased linearly in the initial stage of the SMD simulation. Afterward, at around 9.2 ns, the pulling force reached a peak value of just about 1338.65 pN. At this moment, chlorpyrifos was about to leave the active site of the APH. Time dependence of the distance changes for chlorpyrifos and APH in the two SMD simulations are listed in Figure S5. It can be seen that chlorpyrifos was out of APH through P1 about 9.2 ns. The value of the external force then rapidly decreased to zero and featured minor fluctuations near zero. As seen from Figure 5(b), the pulling force exerted on chlorpyrifos through P2 increased linearly in the initial stage of the SMD simulation until around Table 1. Bottleneck residues of the top ranked tunnels of APH identified by CAVER 3.0 in molecular dynamics trajectory using the probe radius of 10 nm and the clustering threshold of 8.5.  7.7 ns, with a value of 1812.72 pN. The variation of the pulling force reflects the motion of chlorpyrifos in the solvent, proving its full dissociation from the enzyme. Due to the structural character of APH and the molecular size of chlorpyrifos, its main unbinding pathway seemed to lack other possibilities. As discussed above, the main successful unbinding pathway was P2, which is located between blades 1 and 2 in the β-propeller domain (see Figure S7). At the same time, for the purpose of calculating the PMF, the above-determined unbinding pathway, a reaction coordinate, was artificially separated into 12 equal sections. The calculated PMF profiles for P2 and P1 are depicted in Figure 6(a) and (b). The free energy curve reveals some interesting information about the unbinding pathway of chlorpyrifos. With the departure of the inhibitor from the initial equilibrium position, the free energy value rapidly increased, and this is attributed to the interaction of the inhibitor with the active site residues. As is seen from Figure 6(a), the initial position of the inhibitor was in the most stable binding state with the APH. Conversely, in Figure 6(b), with lower distance, the energy curve sharply decreased, and with the departure of the inhibitor from the equilibrium position, the free energy value rapidly increased, which can be attributed to the interaction of the substrate with the active site residues. From 3.0 nm, the energy curve shows equilibrium. It is worth noting that the free energy of unbinding should be calculated in order to obtain a precise, reliable estimate for the relative preference of the ligand unbinding pathways. The free energy barrier of the inhibitor unbinding process for P2 was around 32.74 kcal mol −1 while for P1 it was around 41.00 kcal mol −1 . Thus, the PMF supports P2 as a most favorable chlorpyrifos unbinding pathway.

Unbinding along P2
The above PMF calculation provided a great deal of important information on energy changes during the binding or unbinding process of the substrate, but it is still unknown which atomic interaction results in these energy changes. To explore the atomic essence underlying the energy changes, we carefully investigated the SMD simulation trajectories of the inhibitor's dissociation from APH.
At the beginning, in order for the benzene ring of chlorpyrifos to get through the tunnel, a high force was required to break the hindrance from Arg526. After 900 ps, the Pi cation between Arg526 and chlorpyrifos was broken.
Between 900 and 1000 ps, the force maintained a transient stability. In addition, the middle part of the inhibitor formed van der Waals contacts with Arg526, Phe485, Gly369, Met477, and Trp474, and formed electrostatic interactions with Gly368, Pro370, Ser445, and His556. Between 2000 and 2100 ps, the force maintained a transient stability. The middle part of the inhibitor formed van der Waals contacts with Phe485, Val471, Met477, Thr527, Ala557, and Asp524, and formed electrostatic interactions with Pro370, His556, and Arg526. These strong interactions prevented the inhibitor from releasing quickly (Figure 7(a)).
The force peak at 9220 ps was attributed to the hydrogen bonds of the inhibitor with Gly86 and Asn65 (see Figure 8(a) and (b)), and strong electrostatic interactions with Ile558, Asn559, Ala557, Lys85, Ala87, Glu88, Ser84, Gln89, Asp82, and His90. In addition, Asn65, Asp82, Glu88, Ala557, Ile558, and Asn559 located between blades 1 and 2 in the β-propeller domain were well in line with observations of the crystal structure (Bartlam et al., 2004) (see Figures S6 and S7). From Figure 9(a) and (b), the hydrogen bond between Gly86 and the inhibitor was formed around 4400 ps and broken at 9220 ps, while the hydrogen bond between Asn65 and the inhibitor was formed around 6700 ps and broken at 9220 ps, when the inhibitor was about to leave the APH. Arg526 and Glu88 have been found to be Figure 6. The PMF profile along the unbinding pathway of chlorpyrifos via (a) P2, (b) P1, respectively. As the reverse process, the free energy profile also applies to the binding process.
important residues in stabilizing the ligand binding in the APH pocket (Han et al., 2013;Lu & Chen, 2012). For the rotation, the variation of the between Arg526 and Glu88 was measured during the SMD simulation, as shown in Figure 9(a) and (b). The distance underwent a significant change at 9220 ps. Glu88 had strong electrostatic interactions with the inhibitor and may have induced the side chain with the positive charge of Arg526, to undergo a displacement, partly hindering the inhibitor from leaving. Therefore, Arg526, Glu88,  Asn65, and Gly86 are important residues in the ligand leaving via P2. Feng et al. found that single mutation at position 526 of apAPH can have very different effects on the peptidase and esterase activities with Ac-Leu-pNA and pNPC8 as substrates (Wang et al., 2006). In 2009, Feng et al. also pointed out that in E88A and E88A/R526K mutants, with a broken interdomain salt bridge and a positive charge at position 526, catalytic activities for both a peptidase substrate and an esterase . Our results were consistent with the experimental data.
After chlorpyrifos had gotten rid of the constraints from the residues of the active triad, it could exit smoothly. The fluctuations between 9220 and 10,000 ps were due to the interactions of the inhibitor with the residues in the APH surface. The increase in the force at 9220 ps was because of the breakage of the hydrogen bond formed between Gly86, Asn65, and the inhibitor.

Unbinding along P1
In the first 900 ps, the chlorpyrifos unbinding event along P1 was similar to that along P2. From 910 to 2700 ps, the force decreased slightly, resulting in less blockage by the residues (Phe371, Phe153, Arg526, and Leu115) of interactions with the inhibitor during this period (Figure 10(a)). The force peak formed at 7700 ps, which was attributed to the hydrogen bond of the inhibitor Asp69 and the strong Pi cation interaction with Arg287. In addition, the middle part of the inhibitor formed van der Waals contacts with Phe155, Ser26, and Val118, and electrostatic interactions with Pro370, His556, and Arg526. Ser199, Ser157, Val156, Gly117, Leu27, and Gln28 (Figure 10(b)). From Figure 10(c) and (d), the hydrogen bond between Asp69 and the inhibitor was formed around 5700 ps and broken at 7700 ps, when the inhibitor was about to leave APH via P1. Therefore, Asp69 is important residue for the ligand leaving via P1. The fluctuations between 7700 and 10,000 ps were due to the interactions of the inhibitor with the residues in the APH surface. The increase in the force at 7700 ps was because of the breakage of the hydrogen bond formed between Asp69 and the inhibitor.
The dissociation of chlorpyrifos along P2 and P1 is parallel to its binding orientation in the APH binding  pocket, as shown in Figure 11. Although APH pathways were reliably identified by CAVER 3.0, with estimated relative importance P1 > P2 due to size. Seen form Table S1, the average values of Fmax and Fsum of P2 were lower than that of P1. And so P2 still proved the main tunnel for chlorpyrifos unbinding escape in SMD simulations.
CAVER 3.0 estimated relative importance P1 > P2 for its size (seen from Figure S4 it can seen the size of P1 is a bit larger than that of P2). Rao et al. indicated that the P1 entrance to the active site is via a tunnel in the β-propeller domain, whose diameter is approximately .7 nm (Bartlam et al., 2004). P2 side opening tunnel provides access to the active site about .6 nm (Bartlam et al., 2004). Our results is consistent with the experimental data. SMD can provide the information between the ligand and the residues in the ligand release pathway (Vashisth & Abrams, 2008;Xu et al., 2003;Zhang et al., 2013). P1 entrance to active site is via a tunnel in β-propeller domain, which contains a hydrophobic environment and is particularly rich in phenylalanines (Phe153, Phe155, and Phe163). These residues located near the P1 tunnel and thence can hinder the ligand to release. To sum up, the main tunnel for the chlorpyrifos unbinding escape pathway may be decided to two factors: the size of tunnel and the hinder from the ligand to release. The size of P1 is a bit larger than that of P2 (the diameter: .7-.6 nm; see Figure S6), and so the hinder from the ligand to release will become the main factor for ligand escape. P1 entrance to active site is via a tunnel in β-propeller domain, which contains a hydrophobic environment and hence can hinder the ligand to release. Recently, several mechanisms of POP enzymes were proposed for the access of the ligands (substrates and inhibitors) to the catalytic site across the cavity in the β-propeller and a flexible loop at the interdomain interface, identified in the experimental methods such as apAPH X-ray structure (Bartlam et al., 2004), electron microscopy of a POP enzyme (Tarragó et al., 2009). Some theoretical works (MD, SMD, and US) suggested the interdomain interface being the main tunnel for POP enzymes (Fuxreiter et al., 2005;Kaszuba et al., 2009;Papaleo & Renzetti, 2012;St-Pierre, Karttunen, Mousseau, Rog, & Bunker, 2011). APH is also one of the member of POP enzymes, and so our SMD result that P2 being the main tunnel for APH is consistent with previous work.

Energies of chlorpyrifos binding to the APH via P1 and P2
It is well known that the binding free energy from MM/ GBSA method will provide semi-quantitative estimate for ligand's affinity with enzyme . And so MM/GBSA methods can be used to identify chlorpyrifos unbinding to the APH via P1 and P2 at the force peak. Table 3 lists the binding free energies and their components for two conformations. As shown in Table 3, it should be stressed that the binding free energies (ΔG bind ) of two interaction modes are both the negative values, suggesting that they are energetically favorable. While it is obvious that chlorpyrifos binding to the APH via P2 at the force peak is lower in energy than chlorpyrifos binding to the APH via P1 at the force peak indicating that chlorpyrifos is the highly probable binding conformation and can be easy to be taken off. From the calculations of MM/GBSA binding free energies, we can draw a conclusion that the binding conformation of chlorpyrifos binding to the APH via P2 at the force peak is energetically more favorable than the conformation of chlorpyrifos binding to the APH via P1 at the force peak.

Conclusion
APH from Hyperthermophilic Aeropyrum pernix K1 belongs to the POP family of serine proteases, and catalyzes the N-terminal hydrolysis of N α -acylpeptides to Figure 11. Chlorpyrifos escape from APH along P2 (cyan) and P1 (green). release N α -acylated amino acids. We have applied SMD simulations to identify potential escape routes of chlorpyrifos from hydrolase hydrophobic cavities in the APH-inhibitor complex. Although APH pathways were reliably identified by CAVER 3.0, with an estimated relative importance P1 > P2 for its size, P2 still proved to be the main tunnel for the chlorpyrifos unbinding escape pathway in SMD simulations. Arg526, Glu88, Asn65, and Gly86 are important residues for the ligand leaving via P2. Our study will provided a theoretical clue for APH further study.

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

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work was supported by Major scientific research projects of Jilin Province [20140203025NY], and was worked at the High Performance Computing Center of Jilin University.