Insight into the binding modes and inhibition mechanisms of adamantyl-based 1,3-disubstituted urea inhibitors in the active site of the human soluble epoxide hydrolase

Soluble epoxide hydrolase (sEH) is a promising new target for treating hypertension and inflammation. Considerable efforts have been devoted to develop novel inhibitors. In this study, the binding modes and interaction mechanisms of a series of adamantyl-based 1,3-disubstituted urea inhibitors were investigated by molecular docking, molecular dynamics simulations, binding free energy calculations, and binding energy decomposition analysis. Based on binding affinity, the most favorable binding mode was determined for each inhibitor. The calculation results indicate that the total binding free energy (ΔGTOT, the sum of enthalpy ΔGMM-GB/SA, and entropy −TΔS) presents a good correlation with the experimental inhibitory activity (IC50, r2 = .99). The van der Waals energy contributes most to the total binding free energy (ΔGTOT). A detailed discussion on the interactions between inhibitors and those residues located in the active pocket is made based on hydrogen bond and binding modes analysis. According to binding energy decomposition, the residues Asp333 and Trp334 contribute the most to binding free energy in all systems. Furthermore, Hip523 plays a major role in determining this class of inhibitor-binding orientations. Combined with the results of hydrogen bond analysis and binding free energy, we believe that the conserved hydrogen bonds play a role only in anchoring the inhibitors to the exact site for binding and the number of hydrogen bonds may not directly relate to the binding free energy. The results we obtained will provide valuable information for the design of high potency sEH inhibitors.


Introduction
Soluble epoxide hydrolase (sEH), a member of the α/β-hydrolase fold family of proteins, is the primary enzyme that converts the epoxyeicosatrienoic acids (EETs) to their corresponding dihydroxy eicosatrienoic acids by the addition of water molecule (Newman, Morisseau, & Hammock, 2005;Spector & Norris, 2007;Spector, 2009). EETs exert many biological effects including attenuation of inflammation, vascular smooth muscle cell migration, platelet aggregation, and promotion of fibrinolysis (Behm, Ogbonna, Wu, Burns-Kurtis, & Douglas, 2009;Sun et al., 2002). Therefore, increasing concentrations of EETs by inhibition of sEH may improve the biological functions of EETs. It has been postulated that sEH enzyme can be taken as a potential therapeutic target for cardiovascular diseases (Imig, 2006) and acute inflammation (Schmelzer et al., 2005). In the past decade, great efforts have been devoted to developing the molecular and pharmacological inhibition of sEH in academia, biotechnology laboratories, and pharmaceutical industry (Huang et al., 2012;Imig, 2006;Kim, Morisseau, Watanabe, & Hammock, 2004;Pecic, Deng, Morisseau, Hammock, & Landry, 2012;Rose et al., 2010;Wang, Davis, Jiang, Zhao, & Xu, 2013). New drug design and discovery of highly potent and selective inhibitors of sEH have attracted more attention in many research establishments (Huang et al., 2010). Fortunately, many potent and selective sEH inhibitors including ureas, amides, and carbamates have been synthesized and assessed for their inhibitor activity in vitro and in vivo models (Hwang et al., 2011;Kim et al., 2004;. The primary inhibitors of sEH are epoxides acting as substrate mimics of sEH. However, these compounds are rapidly deactivated by glutathione and glutathione transferases in vivo conditions (Shen, 2010;Spector, Fang, Snyder, & Weintraub, 2004). A significant advance is that ureas, amides, and carbamates were regarded as broad series of transition state analogs of sEH enzyme (Pecic et al., 2012). Urea-based compounds are the most explored sEH inhibitors to date. There are several potent urea-based sEH inhibitors with IC 50 in the low nano-to pico-molar concentrations . Especially, 1,3-disubstituted ureas and related compounds are a class of potent and stable inhibitors of sEH. However, poor solubility, stability, and high melting points of the urea-based inhibitor probably affect their efficacy in vivo (Anandan et al., 2011;Yu et al., 2000). Hence, considerable efforts have been made in improving water solubility and stability of inhibitors as drug candidates. When one adds polar groups to lipophilic inhibitors with the purpose of improving water solubility, the biological activity of the inhibitor usually is reduced. It is a challenge to improve the solubility of those urea-based inhibitors without loss of inhibition potency at the same time. Bruce Hammock developed the so called "Pharmacophore Model" for sEH inhibitor structure-activity relationship (SAR) (Kim et al., 2004, which speeded up the discovery of novel inhibitors (Ingraham, Gless, & Lo, 2011). The "pharmacophore model" consists of three parts, the "primary pharmacophore", the "secondary pharmacophore", and the "tertiary pharmacophore". The "primary pharmacophore" may be a carbamate, amide, or urea, bearing a bulky hydrophobic substituent; the "secondary pharmacophore" should be a polar group, which may form hydrogen bonds to the enzyme within the active site (Kim et al., 2004). Many researches of SAR have indicated that incorporating polar functional groups as the secondary or tertiary pharmacophore into urea-based scaffold is an effective means to improve the water solubility of ureabased inhibitors (Huang et al., 2010). More recently, Huang et al. (2012), successfully optimized and synthesized piperazine-containing 1,3-disubstituted ureas with sub-nanomolar IC 50 and favorable water solubility. What we are interested in is the binding modes and inhibition mechanisms of these novel inhibitors in the active site of sEH.
It is important to understand the binding mode of potent inhibitor in the active site of sEH (Chen, Zhang, Li, & Han, 2012). However, accurate prediction of the binding mode of an inhibitor in the sEH active site is difficult (Li, Chen, Zhao, & Han, 2013). Even for those structurally similar inhibitors, there may still be much difference in the binding mode. For instance, X-ray crystal structures of sEH complexed with four different dialkylurea inhibitors bearing pendant carboxylate "tails" of varying lengths exhibit diverse binding modes, which include two opposite orientations of inhibitor in the active site (Gomez, Morisseau, Hammock, & Christianson, 2006). Moreover, N,N′-diadamantylurea inhibitor possesses high potency against human sEH, which indicates that both sides of the active site tunnel have sufficient volume to accommodate the adamantyl group (Morisseau et al., 2002). Since the bulky adamantyl group can occupy one of two sides of the active site tunnel, it is possible that 1-adamantyl-urea-based inhibitors bind in the active site of human sEH with different orientations. Gomez et al. emphasized that inhibitor-binding orientation in the active site of human sEH cannot be predicted by inhibitor potency (Gomez et al., 2006). Although various inhibitors had been docked into the active site, which helps elucidate the inhibitor SAR (Das, Dhanawat, Kulshrestha, & Shrivastava, 2011;Hwang et al., 2011;Jones, Tsai, Do, Morisseau, & Hammock, 2006;Pecic et al., 2013), there is a limitation of molecular docking due to the accuracy of scoring function.
In order to better examine the principal interactions between adamantane-based 1,3-disubstituted urea inhibitors and key residues of sEH, we calculated and analyzed the binding energy of the reported inhibitors by performing MD simulations. Here, the interactions between the functional groups and residues are described, which may provide a new insight into the structure-affinity relationship.
Understanding the protein-ligand interaction at the atomic level will drive further structure-based drug design processes in the quest for more selective and effective sEH inhibitors.

Ligands preparation
To study the receptor-ligand interaction, five representative adamantyl-based 1,3-disubstituted urea inhibitors were selected from previous experimental results (Huang et al., 2012). In order to compare the distinct inhibitorbinding orientation of two analogues with different methylene units in length, we extracted the two ligands from the crystal structures of enzyme-inhibitor complexes (PDB code: 1ZD2 and 1ZD3) as well (Gomez et al., 2006). All chemical structures of seven inhibitors are shown in Table 1. The labels of all ligands are consistent with relevant publications (Huang et al., 2012;Jones et al., 2006;Kim et al., 2004;Morisseau et al., 2002). All the structures of ligands, except CU2 and CU4, were converted into 3D structure by the GlycoBioChem PRODRG (Schuttelkopf & van Aalten, 2004), and then were fully optimized using Gaussian 09 program at B3LYP/6-31++G(d, p) level (Frisch et al., 2009). Finally, the optimized geometries were used for the subsequent molecular docking.

Construction for the 3D model of receptor-ligand complex
The initial 3D structures of human sEH complexed with CU2 and CU4 were taken from the Protein Data Bank (PDB code: 1ZD2 and 1ZD3). The crystal structure includes a C-terminal domain with epoxide hydrolase activity and an N-terminal domain with phosphatase activity. All crystal water molecules were removed from the PDB file. Our goal is to explore the binding modes and inhibition mechanisms of the inhibitors in C-terminal domain. Hence, only the C-terminal domain is selected for present work and the N-terminal domain was deleted from the PDB file. This would be reasonable because the catalytic activity of two domains is independent of each other (Cronin et al., 2003). The 1ZD3 crystal structure which removed the ligand of CU4 and N-terminal domain was selected as a receptor for subsequent molecular docking.
For docking, the program of AutoDock 4.2 (Huey, Morris, Olson, & Goodsell, 2007) was used for the receptor and ligand docking. The PDBQT files of the receptor and ligands were prepared by AutodockTools 1.5.4. (Sanner, 1999) Meanwhile, all hydrogen atoms, Gasteiger charges, and rotatable bonds were assigned and identified during PDBQT file preparation. For the AutoGrid calculation (Morris et al., 2009), the grid maps employed a grid module with 80 Â 80 Â 80 points of .375 Å spacing. The centroid of the inhibitor was set as the center of the box. For the AutoDock (Morris et al., 2009) calculation, the receptor and ligand docking used Lamarckian genetic algorithm with the following settings of parameters: population size of individuals 150, a maximum number of 25 million energy evaluations, a maximum number of generations of 27,000, a crossover rate of .8, and a mutation rate of .02 were set up. All other docking parameters were set to default. For each docking calculation, 100 independent docking runs were produced. The docked conformations were clustered using a tolerance of 2 Å root-mean-square deviations and evaluated depending on scoring function.
After the docking experiment, 100 independent docking runs were obtained for each inhibitor. The docked conformations of each inhibitor were ranked into clusters according to scoring function. The best conformation was determined based on two criteria, the AutoDock scoring function and both hydrogen bond distances with a cut-off of 2.8 Å. One hydrogen bond is formed between oxygen atom of the urea (or amide) moiety of the inhibitor and the hydroxyl oxygen atom of the side chain of the residue Tyr381 or Tyr465. The other is formed between one of the nitrogen atoms of the urea moiety of the inhibitor and one of the oxygen atoms of the residue Asp333. Concretely, the selection consists of two steps. Firstly, inspection of both hydrogen bond distances for each RMSD cluster. And then, the conformation of the lowest binding energy among several similar docking poses selected from the first step was used as a starting geometry for MD simulation. As both sides of the active site tunnel can accommodate the adamantyl group (Gomez et al., 2006;Morisseau et al., 2002), we obtained two starting geometries with different orientations for each inhibitor at last. According to adamantyl group in the different pockets, two poses of each ligand were divided into two groups (pose-A and pose-B). For pose-A, the adamantyl group lies in a hydrophobic pocket "W334"(constituted by W334, M337, I361, P362, A363, M367, P369, F379, M468, W472, A475, V497, L498, and M502); for pose-B, the adamantyl group locates in an adjacent hydrophobic pocket F265 (including F265, P266, F385, L395, L406, L416, M418, L427, F428, F496, and W524). As a result of structural symmetry of N-N′-adamantyl urea, only one conformation pose was selected. Finally, 11 complex structures were used as starting structure for the following MD simulations. Nine of them were gained from docking results and the other two complexes were obtained from crystal structure.

MD simulations
The initial conformations were acquired from the crystal structure and the best docking results, respectively. Eleven MD simulations, sEH-4c-A, sEH-4c-B, sEH-4d-A, sEH-4d-B, sEH-8d-A, sEH-8d-B, sEH-APAU, sEH-14, sEH-CU2, and sEH-CU4 complexes are performed using the AMBER10 suite using the Amber 03 force field (Case et al., 2008). Amber general force field (Wang, Wolf, Caldwell, Kollman, & Case, 2004) was used for those inhibitors. The missing hydrogen atoms of the proteins were added by using tleap module in the AMBER10 suite. The electrostatic potential of inhibitor was calculated at HF/6-31G ⁄ level using Gaussian 09 package and partial atomic charges were determined by the restrained electrostatic potential (Cornell et al., 1996) fitting technique using Antechamber encoded in AMBER10 suite (Case & Darden, 2008). Protonation states of ionizable residues were determined by H++ (Gordon et al., 2005) at a neutral pH. The final protonation state was as follows: for 1ZD3 crystal structure, HIE237, HIP249, HID263, HIP332, HIE419, HID505, HIE512, HIE517, and HIP523; for 1ZD2 crystal structure, HIE237, HIP249, HID263, HIP332, HIE419, HIE505, HIE512, HIP517, and HIP523, where HID is protonated on the delta nitrogen, HIE is protonated on the epsilon nitrogen, and HIP is doubly protonated. All complexes were solvated in an octahedron periodic box of TIP3P (Jorgensen, Chandrasekhar, Madura, Impey, & Klein, 1983) water molecules with a 12 Å buffer around the complex along each dimension. Na + ions were added to neutralize the net negative charge of the system.
Prior to MD simulation, energy minimization was performed in two consecutive steps, each of which consists of 10,000 cycles (5000 cycles of steepest descent and 5000 cycles of conjugate gradient minimization). In the first step, a harmonic constraint was applied to restraining protein and ligand with a strength of 500 kcal/ (mol Å 2 ). It means that the water and counter-ions can move freely. This step aims to resolve the clashes between the solute and the solvent molecules. The next step, the entire system was optimized without any constraint to eliminate possible remaining bad contacts in the system. After minimization, the system gradually heated from 0 to 300 K in 50 ps at constant volume and then equilibrated at a constant temperature of 300 K and a constant pressure of 1 atm for 500 ps. In MD simulation with 2 fs time step, the Particle-Mesh-Ewald (PME) method (Darden, York, & Pedersen, 1993) was applied to calculating long-range electrostatic interactions and the SHAKE algorithm (Ryckaert, Ciccotti, & Berendsen, 1977) was applied to constraining bonds involving hydrogen atoms. The non-bonded interactions were set to a cut-off distance of 10 Å. Periodic boundary conditions were applied to all dimensions. MD trajectories were produced by PMEMD module and the obtained MD trajectories were analyzed by the PTRAJ module of AMBER. All figures were done with the Chimera program (Pettersen et al., 2004).

MM-GB/SA binding free energy calculations and binding free energy decomposition
The binding free energies between the sEH and the inhibitor were calculated using MM-GB/SA methodology (Case et al., 2005;Feig, Im, & Brooks, 2004;Kollman et al., 2000). To assess the rationality of the sampling strategy of 11 complexes, the RMSD values of backbone atoms of sEH are calculated based on the starting structure. As the flexibility of terminal residues would affect the actual RMSD values, only the residues (N231-A546) were taken into account for the RMSD calculations. As is seen from Figure S1 of the Supplementary Material, all the RMSD curves of complexes under dozens of nanoseconds MD simulations have reached a plateau. In accordance with the RMSD curves, each system remains in equilibrium phase over 20 ns. Hence, it is sufficiently stable and reliable for the subsequent binding free energy calculation and free energy decomposition. One thousand snapshots were extracted from the last 20 ns MD trajectory with 20 ps interval for each system. The water molecules and counter-ions in MD trajectory were stripped by the PTRAJ module of AMBER. For each snapshot, the relative free energy (ΔG bind ) of the inhibitor to the protein was calculated by the follow equation: where G complex ; G receptor ; and G ligand are the free energies of the complex, receptor, and ligand, respectively. Each free energy term in Equation (1) was estimated as a sum of the three terms: the absolute free energy in the gas phase (DG MM ), the solvation free energy (DG sol ), and the change of conformational entropy on ligand term (-T DS), using the equation as follows: The free energy in the gas phase (DG MM ) was further split into electrostatic (DG ele ) and van der Waals (DG vdw ) energies: The solvation free energy DG sol ; can be partitioned into an polar solvation free energy (DG GB ) and a nonpolar solvation free energy (DG SA ) The polar solvation free energy (DG GB ) to the solvation energy was computed with a GB module (Onufriev's GB (Onufriev, Bashford, & Case, 2004), IGB = 2) of the AMBER suite. In our calculations, the dielectric constants were set 1.0 to the interior solute and 80.0 to the exterior solvent, respectively. Atomic radii and charges were in accordance with those used in the MD simulations. The non-polar contribution DG SA was determined according to the Equation (5): where γ, the surface tension proportionality constant, was .0072 kcal/(mol Å À2 ), and β, the offset value, was 0. SASA, determined with the LCPO model (Weiser, Shenkin, & Still, 1999), was estimated by the MSMS algorithm (Sanner, Olson, & Spehner, 1996) with a probe radius of 1.4 Å. The solute entropy contribution (ÀT DS) to the binding free energy was computed through a NMA with the NMODE module in the AMBER10 program (Case et al., 2005). The computations were based on 100 snapshots over the last 20 ns of each trajectory with a 200 ps interval.

Analysis of hydrogen bonds
The hydrogen bonds between both receptors and ligands were calculated using the PTRAJ module of AMBER10 over the last 20 ns of simulations. The cut-off used for hydrogen bonding calculation: bond distance and angle was 3.5 Å and 120°, respectively. The percentage of structures in which the hydrogen bond satisfies the conditions of distance and angle is quantified by percentage of occupation.

Results and discussion
Binding free energy calculation In order to determine the preferential modes of binding of selected inhibitors, the MM-GB/SA method was used to compute binding free energy on a statistically relevant structure ensembles collected along the last 20 ns of the MD simulations. The binding mode with lower binding free energy from the pose-A and pose-B of each inhibitor was determined as the better binding conformation. The calculated binding free energy (ΔG MM-GB/SA ) and detailed energy contribution (van der Walls ΔE vdW , electrostatic interaction ΔE ele , solvation energy ΔG SA + ΔG GB ), and entropy contribution (ÀTΔS) using the single trajectory MM-GB/SA method are summarized in Table 2, respectively. As is seen from Table 2, one can discern the different binding modes easily by comparing with the binding free energy for each compound except compound 8d. For inhibitor 8d in pose-A and pose-B, they are close to each other in total binding free energy (ΔΔG TOT ) containing entropy contribution (ÀTΔS). If the entropy contribution was neglected, the sEH-8d-A mode would be a little more favorable. Both two binding modes may be possible in this system, because we have observed that CU2 and CU4 bind to sEH with opposite orientations. Perhaps, a high level calculation may identify the more favorable binding mode for 8d inhibitor, but it will not be impacted for us to clarify the interaction mechanism. Based on the total binding free energy, we can determine the more favorable binding modes and total binding free energy for each inhibitor. Moreover, when we ranked these inhibitors according to the ΔΔG TOT , we found that the rankings of these ΔΔG TOT remain the same as those of the experimental IC 50 values. The ΔΔG TOT is well correlated with the experimental IC 50 values. The correlation coefficient (r 2 ) reaches up to .99 ( Figure S2 of the Supplementary Material), indicating that theoretical results agree very well with those experimental data. Therefore, the MM-GB/SA method can be used to predict drug potency of this system. In order to better understand which energy term contributes more to binding affinity, the energy components (ΔE ele , ΔE vdW , ΔG SA , ΔG GB , ΔG ele , and ÀTΔS) are compared carefully. Among these energy components, the van der Waals energy (ΔE vdW ) is significantly stronger than the others. Hence, the hydrophobic interactions between ligand and protein are predominant. Additionally, the magnitudes of ΔE vdW are in direct proportion to the size of ligand and adamantyl group which may increase the favorable hydrophobic interaction. The lesser ΔE vdW is the principal reason that CU2 and CU4 present poor potency. The ΔE ele contributes about half of the ΔE vdW to the binding free energy (ΔG MM-GB/SA ). It is worthy to note that the Coulombic interactions are weakly favorable for sEH-4c-B complex formation (ΔE ele , À8.14 Kcal/mol) and are unfavorable for sEH-CU4 complex (ΔE ele , 10.64 Kcal/mol). Since the formation of complex needs to overcome the cost of entropy and desolvation of protein and ligand, polar solvation energy (ΔG ele ) and entropy show an unfavorable contribution for all complexes. Yet, the molecular mechanics energy prominently facilitates complexes over the unbound molecules. Consequently, both electrostatic interactions and van der Waals interactions are crucial to the binding free energy. The entropy change (ÀTΔS) is less compared with the enthalpy change (ΔG MM-GB/SA ). Thus, the binding process that prompts the inhibitors binding to sEH is driven by enthalpy change.

Hydrogen bonds and structure binding mode relationship analyses
In this section, we focus our attention on discussing the interaction characteristics of the different binding modes for each inhibitor in the binding site of sEH. For convenience, a clustering procedure based on the average-linkage algorithm (Jain, Murty, & Flynn, 1999;Kelley, Gardner, & Sutcliffe, 1996) was carried out over the last 20 ns of the trajectory to select representative conformation.
Figure 1(A) shows the overlaps of the crystal structure and the representative structure of sEH-CU2 complex. One can see that a loop region located at the entrance of the binding pocket changes slightly during the MD simulation. This is due to the factor that the loop region is fully exposed in solvent which reduces the interaction with surrounding residues. Therefore, it exhibits high flexibility. However, the flexible regions may become an on-off switch to control the binding pocket in the open and closed conformations. It is advantageous for substrate binding and releasing. In the binding site, the location of the inhibitor CU2 is consistent with that of CU2 in the corresponding crystal structure. Nevertheless, the conformation of CU2 undergoes a slight change during MD simulation. The orientation of urea functional group of inhibitor has approximately 90°r otation, which improves the poor orientation of the hydrogen bond between urea and Asp333 and Tyr465. Additionally, the carboxylate oxygen (O3) of CU2 relative to it in the crystal structure points toward the hydroxyl oxygen atom (OG1) of the polar amino acid Thr358 with a hydrogen-bond distance 2.8 Å. To investigate the hydrogen bonding interaction in detail, hydrogen bonds are calculated along the last 20 ns of MD simulation. Only hydrogen bonds having an occupancy more than 40% during MD simulation were listed in Table 3. It is believed that these hydrogen bonds play crucial roles in the ligand binding. The information for the hydrogen bond with less than 40% occupancy was also listed in Table S1 of the Supplementary Material. From Table 3, we can find that the rotation of urea moiety results in forming three strong hydrogen bonds between CU2 and residues (Asp333 and Tyr465). Therefore, the formed hydrogen bonds indicate that the rotation of urea moiety is in favor of increasing the stability of binding. Three stable hydrogen bonds may help CU2 to bind in active side.
Although there is a subtle difference between the structures of CU2 and CU4, the binding modes of sEH-CU2 (PDB: 1ZD2) and sEH-CU4 (PDB: 1ZD3) are remarkable differences. Figure 1(B) shows the superimposition of the representative structure and the crystal structure of sEH-CU4 complex. It can be observed from the figure that there are no large conformational changes in sEH protein. The conformation of CU4 remains almost the conformation in the crystal structure except for a slight movement toward F265 pocket. The carbonyl oxygen and nitrogen atoms of urea moiety form three stable hydrogen bonds with the side chains of Tyr381 and Asp333, respectively (Table 3). These hydrogen bonds prompt inhibitor to bind on the active site. With respect to the interactions with the surrounding residues, terminal cyclohexane of CU4 makes favorable van der Waals interactions with the indole moiety of Trp334. However, the short aliphatic carbon chain as a linker between urea moiety and the terminal carboxylic group presents a curve because of flexibility, which reduces contact area of van der Waals with surrounding residues. As is shown in Figure 1(B), the short aliphatic carbon chain only contacts with the residues Phe265, Val497, and Leu406. Furthermore, the carboxylic group of CU4 establishes an extremely unstable hydrogen bond with the Hip523 (14.44% occupancy, Table S1). From the calculated components of binding free energy, the van der Waals interactions (ΔE vdw ) are the primary factor that promotes the CU4 binding to the active pocket, not the electrostatic interaction (ΔE ele ) due to the contribution of 10.64 kcal/mol. Two factors may be adduced to explain the relatively poor inhibitory potency of CU2 and CU4. One is that the short aliphatic carbon chains of CU2 and CU4 reduce van der Waals contacts with those residues in the active site. Another is that the free acid functionality at the end of short aliphatic chains negatively impacts inhibitory potency. Despite the poor inhibitory potency of CU2 and CU4, they can bind to the active site of sEH indicating that both sides of hydrophobic pocket have sufficient cavity volumes for accommodating bulky hydrophobic pharmacophore.
Morisseau et al. (Morisseau et al., 2002) investigated the potency of the N,N′-diadamantylurea inhibitor (labeled 14) which exhibits high potency against human sEH (IC 50 = 6.4 nM). It may be speculated that both binding pockets of sEH could accommodate the bulky adamantyl group. To verify the speculation, dynamical situation of the N,N′-diadamantylurea inhibitor in the active site has been evaluated by performed MD simulation. The superimposition of the initial structure and representative structure of sEH-14 complex is provided in Figure 2. The superimposed structure reveals that α helix domain consisting of L498-E507 has a greater movement compared to the initial structure. It causes the entrance of the W334 pocket turned off. With careful inspection of the residues around adamantyl group in the W334 pocket, it was found that residue Leu498 from the crystal structure is close to adamantyl group. In order to obtain the optimal contact distance from the side chain of Leu498, α helix domain shifts slightly away from adamantyl group during the energy minimization phase. Adamantyl group leads to receptor conformational changes which make the entropy increasing, but this can be compensated by increasing van der Waals interactions. Furthermore, a steric clash is found between the CE1 atom of residue Tyr381 in the crystal structure and the closest C atom of the adamantyl group in F265 pocket (clash distance: 2.7 Å). To avoid a steric clash, the side chain of Tyr381 moves toward the other side resulting in loss of hydrogen bonding interactions with carbonyl oxygen. Although bulky adamantyl group affects conformations of the active site, the binding free energy indicates that extensive hydrophobic contacts compensate the adverse impact caused by conformational changes.
Many researches have shown that AUPU (AR9281) is a potent and selective inhibitor of sEH in many models (Anandan et al., 2011;Chen, Whitcomb et al., 2012;Chen, Zhang et al., 2012;Imig, Carpenter, & Shaw, 2009). Figure 3 shows two quite different binding modes of AUPU. For better illustrating two binding modes, the conformations with distinct binding modes were superposed. As can be seen from the overall superposition, two loop regions located at the entrances of hydrophobic pockets of sEH-AUPU-A complex and sEH-AUPU-B complex (shown in intense color) are exposed in solvent and are subject to the influence of solvent molecules. So they have high flexibility. Binding modes of sEH-AUPU-A-A and sEH-AUPU-A-B appear to be obtained by mutually rotating 180°around the carbonyl group (C=O) of AUPU. Since the adamantyl group and the piperidine ring locate at the same position of W334 pocket, both adamantyl group and piperidine ring make hydrophobic contacts with Trp334. Electrostatic interactions may exist between the terminal carbonyl oxygen of AUPU and Hip523 in F265 pocket. According to the hydrogen bond analysis (Table 3) the urea group forms stable hydrogen bonds with side chains of Asp333 and Tyr381 in both binding modes. Additionally, the carbonyl oxygen of the urea group can also form stable hydrogen bond with Try465 in sEH-AUPU-A-A complex, while the corresponding carbonyl oxygen in sEH-AUPU-A-B complex loses hydrogen bond with Try465 due to the residue shifts away. Figure 4 presents the superposition of representative conformations of sEH-4c-A and sEH-4c-B complexes. It is seen from Figure 4 that receptor structures display similar geometrical feature during MD simulation. However, the primary pharmacophore (urea group) exhibits noticeable differences in sEH-4c-B binding modes. For sEH-4c-A complex, the urea group engages six hydrogen bond interactions with Asp333, Tyr381, and Try465 (Table 3). Moreover, the adamantyl group forms a π-σ interaction with the indole ring of Trp334 and the benzene ring formed a π-π stacking interaction with the imidazole ring of the His523. In the case of sEH-4c-B binding mode, the urea group loses almost all hydrogen bond interactions with residues (Table S1) because the situation changed. It seems that the length of the W334 pocket is insufficient for accommodating the 4c inhibitor. A probable reason is that tertiary (piperazine ring) and second pharmacophore (benzoyl) are rigid structures, Figure 2. The superimposition of the representative structure (spring green) and initial structure (orange) of sEH-14 complex. which warp the structures hard to adapt the pockets or point to the entrance of pocket. Thus, it forces the urea group to move toward to F265 pocket resulting in losing hydrogen bond interaction. In addition, a gap is formed between the urea group and Asp333. By analyzing the solvent molecules around the urea group during the MD simulation, one water molecule at least forms hydrogen bonds with N atoms of urea group and OD1 or OD2 of Asp333. This reduces the electrostatic interaction between urea group and aspartic acid, which implies that it is an unfavorable binding mode. In accordance with these facts, we can speculate that the sEH-4c-A conformation is a more favorable conformation. It coincides with total binding free energy (ΔΔG TOT ). Figure 5(A) shows the superposition of sEH-4c-A, sEH-4d-A, and sEH-8d-A complexes. As for the complex sEH-4d-A, there is no obvious distinction relative to sEH-4c-A complex. The only difference between 4c and 4d is the substitutional group at the end of the inhibitors. So, sEH-4c-A and sEH-4d-A complexes remain in the same binding pattern. A similar result is observed: the urea group of sEH-4d-A complex engages in the hydrogen bond formation with the active site residues Asp333, Tyr381, and Tyr465; the benzene ring formed a π-π stacking interaction (edge-to-face) with the imidazole ring of the His523 residue too. Compared with sEH-4c-A complex, the terminal methyl of 4d compound forms π-σ interaction with phenyl of Phe496 ( Figure S3 of the Supplementary Material); the terminal tert-butyl group of 4c compound has extensive hydrophobic contacts with surrounding residues relative to the sulfonyl group of 4d compound, which may further contribute to the inhibitor's good binding affinity. Although sulfonyl oxygen of 4c inhibitor is near the hydroxyl oxygen of  Ser413, hydrogen bond analysis indicates that only 7.42% hydrogen bonding occupancy between them (Table S1). It is possible that solvent water molecules enter the hydrophobic pocket and involve hydrogen bonding interaction with sulfonyl oxygen disturbing the binding between sulfonyl group and residues. In addition, the secondary and tertiary polarity pharmacophores neither form hydrogen bonds with the surrounding residues nor form hydrogen bonds with solvent molecules that enter into the hydrophobic pocket. The probable reason is that both sides of the polar functional groups of ligand are linked by hydrophobic groups, which prevent solvent molecules from interacting with the polar functional groups of ligand. In sEH-8d-A binding mode, the structure of the receptor is analogous with the conformations of sEH-4c-A and sEH-4d-A except the loops of entrance of the hydrophobic pocket ( Figure 5(A)). Compared with sEH-4c-A and sEH-4d-A complex, the adamantane moiety and amide moiety of 8d inhibitor have similar binding conformation and interaction mechanism with surrounding residues as described above. The major differences concentrate in F265 pocket. Interestingly, since the right side of 8d is a more flexible phenyloxy alkyl chain, the orientations of the tertiary piperazine pharmacophore and terminal methylsulfonyl substituent have largely changed after MD simulation. The terminal methylsulfonyl substituent is pointed to the mouth of the hydrophobic cavity F265 in the starting structure, while the methylsulfonyl substituent inserts deeply into the bottom part of the F265 cavity during MD simulation. In this direction, it may increase the electrostatic interactions between polar pharmacophore piperazine and polar residues (Arg401, Thr402, and Ser405) and hydrophobic interactions between methyl of the methylsulfonyl substituent and non-polar residues (Leu395, Leu406, Leu427, and Phe428) ( Figure S4 of the Supplementary material). Figure 5(B) shows the superimposition of sEH-4c-B, sEH-4d-B, and sEH-8d-B representative snapshots. It can be observed that there are obvious differences among the three binding conformations of the inhibitors. Similarly, sEH-4d-B binding mode ( Figure S3 of the Supplementary Material) could not form a stable hydrogen bond with Asp333 as sEH-4c-B binding mode due to a water molecule locating at the position between urea group and Asp333. In addition, the benzene ring as a linker between the primary pharmacophore and the secondary pharmacophore flips 180 degrees compared with sEH-4c-B. This action leads to the benzene ring forming in π-π interaction with the indole ring of Trp334 and the methylsulfonyl inserting deeply into the bottom of W334 pocket. In the case of sEH-8d-B binding mode, the amide moiety establishes hydrogen bond interactions with Asp333 and Tyr465. An additional π-σ interaction between adamantyl group and benzene ring of Tyr381 is found in sEH-8d-B binding mode ( Figure S4 of the Supplementary material). Unlike sEH-4c-B and sEH-4d-B, the benzene ring of 8d does not remain π-π interaction with the indole ring of Trp334 anymore, but forms Tshaped contact with the indole ring of Trp334. Furthermore, the terminal methyl is orientated towards the hydrophobic cavity entrance and makes a tight hydrophobic contact with residues of the cavity wall (Ile361, Pro362, Ala363, Met367, and Met502).
From the identification of the most energetic favored binding mode of piperazine-containing 1,3-disubstituted urea inhibitors, we obtain many significant characteristics. All binding modes show that the phenyl rings as a linker between primary pharmacophore and secondary pharmacophore are directional. The phenyl rings keep horizontal and form π-π interaction with the imidazole ring of Hip523 in F265 cavity, while keep perpendicular in W334 pocket. Although few solvent water molecules run into hydrophobic pocket, there is no hydrogen bonding interaction between water molecule and secondary pharmacophore or tertiary pharmacophore. Hence, those inhibitors could accept polar groups without reduced inhibitory potency. However, as polar groups are introduced, it may be beneficial for increasing the water solubility and metabolic stability of inhibitors.

Decomposition of binding energy on per-residue
To further examine the residues contribution to the ligand binding, the MM-GB/SA binding energy (the binding enthalpy) decomposition approach was employed to decompose ΔG MM-GB/SA . Based on per-residue binding enthalpy decomposition, we can determine some key residues that contribute to ligand binding. Figure 6 shows the contribution of individual residue to binding enthalpy. For comparison, Figure 6 was lined into two rows according to binding modes. Those pose-A binding modes were put on the left, while pose-B binding modes were put on the right. Figure 6 reveals that the interaction between the inhibitors and sEH mainly locates on two conserved hydrophobic regions (F265 pocket and W334 pocket) and the active site (including Asp333, Tyr381, Tyr465, and His523). From Figure 6, one can see that the major favorable energy contribution derives predominantly from the residues Asp333 and Trp334 of sEH for all systems, while Hip332 contributes positive energy to the binding free energy in all systems (except sEH-CU2 and sEH-CU4) indicating that it is unfavorable for ligand binding. Asp333 contributes most to the binding free energy in all systems with the exception of sEH-4c-B due to loss of the hydrogen bonding interaction with Asp333. The Trp334 has a significant van der Waals contribution among the complexes because of forming π-π interaction with phenyl of the inhibitor or π-σ interaction with adamantyl moiety of the inhibitor. The most interesting find of all is that the residue Hip523 presents a different mechanism of action in pose-A systems and pose-B systems. For pose-A systems, histidine acid 523 contributes negative binding free energy, which plays a role that encourages the inhibitor binding to the receptor. Conversely, the histidine acid 523 in pose-B system always contributes positive binding free energy, which is unfavorable for inhibitor binding to the receptor. The difference may have something to do with the adamantyl group. In the following, we listed the contribution of many key residues to the binding free energy (Table 4). It suggests that the adamantyl group in F265 hydrophobic pocket shows electrostatic repulsive interaction with Hip523. In order to gain better insights into the binding mechanism, we listed van der Waals contribution (ΔΔG vdW ) and the polar contribution ΔΔG ele decomposition of many key residues in Table 4. (ΔΔG ele is the sum of the electrostatic solvation free energy ΔG GB and MM electrostatic energy ΔE ele .) As is expected, only charged residues have a remarkable electrostatic contribution to the binding. It can be found that unfavorable contribution of the Hip332 residue mainly comes from the electrostatic interaction. And then, both van der Waals energy (ΔΔG vdW ) and the electrostatic energy (ΔΔG ele ) of the Asp333 residue contribute negative values in most systems. The non-polar residue Trp334 has significant ΔΔG vdW contribution. Surprisingly, the polar contribution of polar residues Tyr381 and Tyr465, which may form hydrogen bonds with urea moiety or amide moiety of inhibitor, appears to have an unexpected negative effect on inhibitors' binding. The favorable contribution of polar residues Tyr381 and Tyr465 to the binding energy chiefly comes from van der Waals energy. The reason why the values of ΔΔG ele are positive is that the favorable electrostatic energy (ΔE ele ) is offset by the unfavorable electrostatic solvation free energy (ΔG GB ). Thus, the conserved hydrogen bonds only play a role in anchoring inhibitors to the exact site for binding with sEH. In particularly, the hydrogen bonds formed with the Asp333 play a critical role. Combined with the results of hydrogen bonds analysis, we can draw such a conclusion that the number of hydrogen bonds formed between the receptor and the ligand is not directly related to the corresponding binding free energy; however, the binding free energy is related to the strength of hydrogen bond. Since the oxygen of Asp333 and two N atoms of urea group may form tri-dented hydrogen bond, which is stronger than those common hydrogen bonds, the contribution of electrostatic energy of the residue Asp333 is larger than that of the residues Tyr381 and Tyr465.
In addition, the electrostatic contribution of Hip523 in all systems (except sEH-CU4 and sEH-4d-A) are positive, even though the imidazole ring of Hip523 forms π interaction with the phenyl of 4c-A, 4d-A and 8d-A. In view of the complexity of the issue of π stacking, we simply compared the relative position between the imidazole ring of Hip523 and the phenyl ring of the inhibitors 4c-A, 4d-A, and 8d-A ( Figure S5 of the Supplementary material). The Figure S5 shows that all the conformations of π stacking formed between the imidazole ring of Hip523 and substituted benzene ring of 4c, 4d, and 8d are parallel-displaced conformations. However, high-level ab initio calculations indicates that the perpendicular T conformation may be more stable than parallel-displaced conformation and the energy of an aromatic bond depends on their mutual orientation and distance between aromatic planes (Karthikeyan & Nagase, 2012;Zhuravlev, Shchegolev, Savvateeva-Popova, & Popov, 2009). In contrast, the phenyl in the conformation of sEH-4d-A system has larger movement comparing to the other two. This may be the reason that the ΔΔG ele of Hip523 of sEH-4d-A is lower than the others. Moreover, Karthikeyan (Karthikeyan & Nagase, 2012) reported that the dispersion energy is dominant for the benzene-imidazole and benzene-indole dimers in the total interaction energy. Our results agree with their findings. The repulsive electrostatic energy can be offset by the van der Waals interaction. In conclusion, the Hip523 residue would be the most important residue in determining the orientation of the adamantyl-based inhibitors.
In the end of this section, we focus on comparing the binding affinity in the two hydrophobic pockets (namely, F265 pocket and W334 pocket). The contribution of those residues located in P265 pocket and W334 pocket was summed up, respectively, and was listed in Table S2. It can be observed that the sum of hydrophobic interactions in the W334 pocket is always larger than those in the F265 pocket, regardless of the left part or the right part of the urea/amide group. The length of the inhibitor may influence the binding affinity. Those inhibitors with long hydrophobic chain have more van der Waals contacts with surrounding residues. Furthermore, if the contribution of those polar residues (such as Asp333, Tyr381, Tyr465, and Hip523) in the catalytic sites were not taken into account, the pose-B binding mode is more favorable than pose-A binding mode. From a different perspective, those polar residues play a crucial role in determining the orientation of the ligand, especially, Asp333 and Hip523.

Conclusions
In the present study, molecular docking, MD simulations were performed to investigate binding modes and the interaction mechanisms of a series of highly potent sEHinhibitors. MM-GB/SA free energy calculation was employed to determine the favorable binding mode. According to the ΔΔG TOT , those favorable binding modes are well correlated with the corresponding IC 50 values (r 2 = .99). Decomposition of the binding energy shows that the van der Waals energy is the most predominant factor to derive binding.
Hydrogen bond analysis shows that many important hydrogen bonds already observed in the crystal structure were conserved through the MD simulation. Residues of Asp333, Tyr381, and Tyr465 play a role in stabilizing the inhibitors' association by hydrogen bonding interaction. The interaction features for each inhibitor with adverse binding orientation in the binding pocket are discussed in detail. Some of important conclusions are summarized. Such as: (1) CU2 and CU4 present a poor potency which is due in part to the short aliphatic carbon chain; (2) for selected piperazine-containing 1,3-disubstituted urea inhibitors, the phenyl endows with a definitive direction in the two pocket; (3) there is no hydrogen bond formed between the solvent water molecules and polar groups.
Decomposition of binding free energy indicates that the Asp333 and Trp334 are key residues. Asp333 contributes the most binding free energy in all systems with the exception of sEH-4c-B and Trp334 has a compact hydrophobic contact with inhibitors. The Hip523 shows diverse behavior in different binding modes. For pose-A systems, Hip523 contributes negative binding free energy, while it contributes positive free energy in pose-B binding modes. Combined with the results of hydrogen bond analysis and binding free energy, we believe the conserved hydrogen bonds play a role only in anchoring the inhibitors to the exact site for the binding and the number of hydrogen bonds may not directly relate to the binding free energy.