Understanding the molecular basis of agonist/antagonist mechanism of human mu opioid receptor through gaussian accelerated molecular dynamics method

The most powerful analgesic and addictive properties of opiate alkaloids are mediated by the μ opioid receptor (MOR). The MOR has been extensively investigated as a drug target in the twentieth century, with numerous compounds of varying efficacy being identified. We employed molecular dynamics and Gaussian accelerated molecular dynamics techniques to identify the binding mechanisms of MORs to BU72 (agonist) and β-funaltrexamine (antagonist). Our approach theoretically suggests that the 34 residues (Lys209–Phe221 and Ile301–Cys321) of the MORs were the key regions enabling the two compounds to bind to the active site of the MORs. When the MORs were in the holo form, the key region was in the open conformation. When the MORs were in the apo form, the key region was in the closed conformation. The key region might be responsible for the selectivity of new MOR agonists and antagonists.

Moreover, identifying the ligand-free states of GPCRs can facilitate the development of more selective drugs capable of modulating a specific signaling pathway, thereby minimizing undesirable side effects and improving therapeutic efficacy 10,11 . The X-ray structure of the MOR has been determined in the active state, in which the MOR is bound to the morphinan agonist BU72 12 . Currently, X-ray studies have revealed the inactive structure of the MOR 13 . Furthermore, the X-ray structures of the inactive/active states of the MOR have been obtained 12,13 ; however, because of the lack of experimental conformation of the MOR, many problems remain unresolved.
Studying the binding mechanisms of agonists and antagonists to GPCRs is difficult because long-time scale all-atom dynamics simulations are necessary for sampling conformational states of GPCRs 14,15 . The application of all-atom molecular dynamics (MD) simulations for studying conformational ensembles obtained from a single, long-time-scale conventional molecular dynamics (cMD) simulation is still limited; this limitation is due to the possible energy barriers between various ligand-free states. Thus, an enhanced sampling technique is required for this task. Enhanced sampling techniques have been applied to predict the structural dynamics of GPCRs [16][17][18][19][20][21] . Recent reports show that enhanced sampling techniques have been successfully applied for evaluating binding mechanisms and structural dynamics 17 , including the metadynamics method 22 , adaptive biasing force 23 method and coarse-grained conformational sampling, cMD 14 , and accelerated molecular dynamics (aMD) or Gaussian accelerated molecular dynamics (GaMD) 18 . These enhanced sampling studies provide crucial insights into binding mechanisms and structural dynamics. The disadvantage of enhanced sampling techniques is the requirement of predefined parameters (i.e., root-mean-square distance and protein structures). However, the enhanced sampling method of aMD (or GaMD) avoids such a requirement. In the aMD method, a boost potential is added into the potential energy surface; thus, the energy barriers are effectively decreased, accelerating transitions between the low-energy states 18,24,25 . The aMD method has also been successfully applied to biological system simulations, and aMD simulations conducted on the time scale of hundreds of nanoseconds can approach cMD simulations conducted on the millisecond timescale [26][27][28][29] .
The drawback of the aMD method is the large energetic noise occurring during reweighting 30 . In aMD simulations, the applied boost potential is typically on the order of tens to hundreds of kilocalories per mole (kcal/mol), which is much higher than that of other enhanced sampling methods that utilize protein structures or reaction coordinates. Accurately reweighting aMD simulations is difficult, particularly for large protein molecules 31 . Miao et al. provided a solution (i.e., GaMD) for improving the aMD method. In the GaMD method, the boost potential follows a near-Gaussian distribution, and cumulant expansion to the second order provides improved reweighting of aMD simulations 32 . The reweighted free energy profiles of GaMD are in good agreement with those of the long-time-scale cMD simulations 33 . However, GaMD has the limitation that it cannot evaluate proteins with less than approximately 35 amino acid residues 33 .
In this study, we applied the GaMD method to simulate the binding mechanisms of agonists and antagonists to a MOR and observed the structural dynamics of the MOR.

Results and Discussion
Free energy calculation (PMF) of complex MORs by using GaMD simulations. Free energy (PMF) profiles of complex systems were explored using GaMD simulations of MOR distance values, and the profiles are illustrated in Fig. 1. Snapshots of MORs with agonist (BU72) and antagonist (β-funaltrexamine) ligands are presented in Figures S1 and S2. Our PMF calculations indicate the presence of two energy barriers (major barrier: at RCs of 4-12 Å; minor barrier: at RCs of 18-23 Å) in the five independent 1000-ns GaMD simulations. For the MOR with agonist (BU72), the major energy barrier was 7.19 ± 1.22 kcal/mol and the minor barrier was 2.89 ± 0.68 kcal/mol. For the MOR with antagonist (β-funaltrexamine), the major energy barrier was 6.46 ± 1.06 kcal/mol and the minor barrier was 2.07 ± 0.45 kcal/mol. Moreover the energy barriers were subjected to RMSF calculations, and the snapshots of RCs (3 and 18 Å) were subjected to functionally key residue analysis.
Functionally key residues. Identification of functionally key residues can provide a clear insight into the structural aspects of MORs. In this study, a structure-based approach was applied to identify functionally key residues. According to the snapshots of the RC (18 Å) and the X-ray structures of the MORs, the key residues and pharmacophore regions were analyzed using the Ligandscout program. About the snapshots of the RC (18 Å), the residues (probability that more than half) were selected for the binding mode analysis. Our results are presented in Table 1 Fig. 2.
MORs at an RC of 28 Å. Identifying the apo forms of MORs can provide a clear three-dimensional structure of MORs for designing drugs. Figure S3 shows two apo forms of MORs (at an RC of 28 Å). The RMSD between the two MOR apo forms was 2.05 Å. Figures S3(C) and S3(D) present a comparison of two X-ray structures of MORs with the apo forms of MORs.
Electrostatic and van der Waals binding interactions (major barrier: at RCs of 4-12 Å; minor barrier: at RCs of 18-23 Å). Table 2 shows the electrostatic/van der Waals binding interactions between
Binding mechanism of β-funaltrexamine (antagonist) to MORs. We also observed in our PMF profiles that β-funaltrexamine must overcome the two energy barriers (major barrier: at RCs of 4-14 Å; minor barrier: RC at 18-23 Å) to bind with the binding pocket (     Table 2 showed that the binding interactions were quickly decayed within the RCs of 4-7 Å. Our results showed that Met65-Phe87, Leu116-Ser145, Pro201-Asn230, and Ile301-Cys330 residues might play important roles in relaxing Mors and making BU72 easy to overcome the major energy barrier.

Comparing specific changes of the snapshots of RC at 28 Å and X-ray MORs (RC at 3 Å).
For the binding modes of MORs with BU72 and β-funaltrexamine, our results were presented in Figure S5. At RCs of 18-23 Å, our predicted binding mechanisms showed that no residues interacted with BU72 and β-funaltrexamine, and nine residues (Val126, Asn127, Tyr128, Leu129, Met130, Gly131, Thr132, Trp133, and Pro134) exhibited obvious fluctuations and enabled the two compounds to bind to MORs. At RCs of 4-11 Å, our predicted binding mechanisms revealed that 34 residues (Lys209-Phe221 and Ile301-Cys321) exhibited obvious fluctuations and enabled the two compounds to bind with MORs. Figure S4 illustrates the side and top views of the 34 residues (Lys209-Phe221 and Ile301-Cys321) among the four MORs structures (snapshots at an RC of 28 Å: MOR with BU72 and MOR with β-funaltrexamine; X-ray MOR at 3 Å: MOR with BU72 and MOR with β-funaltrexamine). Our results indicated that the 34 residues (Lys209-Phe221 and Ile301-Cys321) were the key regions enabling the two compounds to bind to the active site of the MORs. Our results indicated that the 34 residues (Lys209-Phe221 and Ile301-Cys321) were the key regions enabling the two compounds to bind to the active site of the MORs. When the MORs were in the holo form, the key region was in the open conformation ( Figure S4, red part) and the BPSAs declined sharply ( Figure S6). When the MORs were in the apo form, the key region was in the closed conformation ( Figure S4, green part) and the BPSAs declined sharply ( Figure S6).

Comparing specific changes of the snapshots of RC at 4 and 18 Å. The snapshots of RC at 4 and 18 Å
were performed for the BPSA analyses by using Hollow and UCSF chimera software 34,35  Comparing the alternative models describing the transition between active and inactive states in GPCRs. Although Prof. Michel Bouvier proposes the hypothesis of the alternative models describing the transition between active and inactive states in GPCRs ( Figure S7) 8 , there is no experimental structural evidence of the ligand-free state in the mu opioid receptor. But, there is a few evidence for ligand-free state of GPCR, such as β1 adrenergic receptor. For the β1 adrenergic receptor, the experimental data (Table S1 and Figure S8) reported by Dr. Huang et al. show that the ligand-free state is similar with active state, and the data also support the two-state model illustrated in Figure S7(B) 36 . Comparing the B1AR ( Figure S8) and MOR ( Figures S3 and  S4), the obvious conformational changes occur in TM1 and TM6, respectively. Dr. Miao et al. used the GaMD to study the activation of M2 muscarinic GPCR 37 . The GaMD method may be suitable for studying the activation of GPCRs with agonist and antagonist. Thus, we used the GaMD to study the binding mechanism of MORs with agonist and antagonist. Finally, our GaMD simulation results tended to Figure S7(A).

Conclusions
In this study, we used GaMD simulations and X-ray structures (agonist and antagonist ligands bound to MORs) to identify the binding mechanisms of MORs to BU72 and β-funaltrexamine. From the X-ray structures, the RCs were defined as the distance between the CM of the compounds and the CM of their binding pockets. Subsequently, we applied the GaMD enhanced sampling method and performed RMSF and PMF calculations to predict the binding mechanisms of the two compounds to MORs. Our PMF calculations indicate the presence of two energy barriers (major barrier: at RCs of 4-14 Å; minor barrier: at RCs of 18-23 Å) in 1000-ns GaMD simulations. For the MOR with agonist (BU72), the major energy barrier was 6.43 kcal/mol and the minor barrier was 1.14 kcal/mol. For the MOR with antagonist (β-funaltrexamine), the major energy barrier was 5.87 kcal/mol and the minor barrier was 1.19 kcal/mol. According to our RMSF profiles, the 34 MOR residues (Lys209-Phe221 and Ile301-Cys321) were the key regions enabling the two compounds to bind to the active site of the MORs. Our results indicated that the 34 residues (Lys209-Phe221 and Ile301-Cys321) were the key regions enabling the two compounds to bind to the active site of the MORs. When the MORs were in the holo form, the key region was in the open conformation ( Figure S5, red part) and the BPSAs were increased ( Figure S7). When the MORs were in the apo form, the key region was in the closed conformation ( Figure S5, green part) and the BPSAs were decreased ( Figure S7). The key region might be responsible for the selectivity of new MOR drugs.

Method
Gaussian accelerated molecular dynamics (GaMD). GaMD is an enhanced conformational sampling method of biomolecules by adding a harmonic boost potential to smoothen the system potential energy surface 32 . when the system potential (V) is lower than a referenced energy (E), a harmonic boost potential (ΔV) is added as: where K is a harmonic force constant. The modified system potential (V*) is given by: 2 IF the system potential (V) is great than a referenced energy (E), a harmonic boost potential (ΔV) is equal to zero.
Smoothening the potential energy surface for overcoming intermedia energy barriers, the boost potential is to satisfy the following step. There are two potential energy values V1 and V2. If V1 < V2, the biased V1* < V2*. By replacing V* with eq. (2), the relationship will be: Step (1) If V1 < V2, the potential difference on the smoothened energy surface should be smaller than that of the original energy surface. By replacing V* with eq. (2), the relationship will be: Step (2) Combing the eq. (3), eq. (4) and the relationship (Vmin ≦ V1 < V2 ≦ Vmax), we can derive: Step (3) Where Vmin and Vmax are the minimum and maximum potential energies. By eq. (5), we can obtain k0 is the magnitude of the applied boost potential.
Step (4) The standard deviation of ΔV must be small enough to ensure accurate reweighting 38 .
where the Vave and σV parameters are the average and standard deviation of the potential energies, and σΔV is the standard deviation of ΔV with σ0 as a user-specified upper limitation for accurate reweighting potential energies. Here the standard deviations of the total potential and dihedral potential boosts are equal to 10 kcal/ mol in our simulations.
SCIenTIfIC RepORTS | 7: 7828 | DOI:10.1038/s41598-017-08224-2 Step (5) Extending the step (2). If E = Vmax, we can drive the eq. (5) and obtain According to eq. (21) and eq. (19), the K0 can be defined as: Step (6) Extending the step (2). If E = Vmin + 1/k, we can drive the eq. (8) and obtain Step (7) GaMD can provide the total potential boost, dihedral potential boost, and the dual potential boost to accelerate the molecular simulations. The boost potential (ΔV) is given as: where K0 is the magnitude of the applied boost potential, Vmin and Vmax are the system minimum and maximum potential energies. The initial K0 is equal to 1.0, and the Vamx and Vmin will be obtained form our cMD simulations. To characterize the extent to which ΔV follows Gaussian distribution, its distribution anharmonicity 32 . GaMD method has been applied in the alanine dipeptide, chignolin and lysozyme simulations 32 .

Reweighted free energy calculations for GaMD simulations (Gaussian Approximation).
The probability distribution of the selected reaction coordinates A(r) is defined as P*(A), where r can be distance, angle, RMSD, etc. 38 . According to the GaMD boost energies of each reaction coordinate, P*(A) can be reweighted and defined as where M is the number of bins, β is equal to KBT, β∆ e V r j ( ) is the ensemble-average factor of the jth bin. For reducing the energetic noise, the ensemble-average factor can be defined as: After driving the eq. (14), the first three cumulants can be defined as: The reweighted free energies can be calculated by GaMD simulation of MORs. Firstly, we modified the inactive MOR pdb file, and we used pymol software to break the covalent bond of the antagonist (β-funaltrexamine) with Lys233 residues. Secondly, we generated our initial models (inactive, PDB ID: 4DKL/our modified pdb file; active, PDB ID: 5C1M) by using the CHARMM-gui server 39 . The initial MOR structures were generated and then inserted into solvent molecules. The solvent molecules contained a POPC lipid bilayer with 20% cholesterol, TIP3 water, and 0.15 M NaCl molecules 40,41 . The size of the MOR system was approximately 11.00 × 11.00 × 14.00 nm 3 . The initial MOR structures were then simulated with the AMBER 14 package by using the AMBER FF14 all-hydrogen amino acid, AMBER lipid 14, and AMBER GAFF force field parameters. The AMBER GAFF partial atomic charges are often based on the RESP fitting procedure of the electrostatic potential obtained at the HF/6-31 G(d) level of theory. The geometries of a morphine agonist (BU72) and antagonist (β-funaltrexamine) were fully optimized, and their electrostatic potentials were obtained using a single-point calculation. Both operations were performed at the HF level with the 6-31 G(d,p) basis set by using the GAMESS US program 42 . All cMD simulations were performed in the isothermal-isobaric ensemble with a simulation temperature of 310 K, unless stated otherwise, by using a Verlet integrator with an integration time step of 0.002 ps and SHAKE constraints 43 of all covalent bonds involving hydrogen atoms. In the electrostatic interactions, atom-based truncation was performed using the PME method 44 , and the switch van der Waals function was used with a 2.00 nm cutoff for atom-pair lists. The complex structure was minimized for 100,000 conjugate gradient steps and was then subjected to a 100-ns isothermal, constant-volume MD simulation and five independent 1000-ns GaMD simulations. Figure 5A shows the initial structure of the active MOR with UB72.