To probe the activation mechanism of the Delta opioid receptor by an agonist ADL5859 started from inactive conformation using molecular dynamic simulations

Abstract The δ-opioid receptor (DOR) is a critical pharmaceutical target for pain management. Although the high-resolution crystal structures of the DOR with both agonist and antagonist have recently been solved, the activation mechanism remains to be elusive. In this study, a DOR agonist ADL5859 was docked to the inactive DOR and multiple microsecond molecular dynamic (MD) simulations were conducted to probe the activation mechanism. While the receptor with the crystal ligand (i.e. antagonist naltrindole) maintained the inactive conformation in all three independent simulations, the receptor with ADL5859 was adopting toward the active conformation in three out of six independent simulations. Major conformational differences were located on transmembrane (TM) 5 and 6, as well as intracellular loop 3. Compared to naltrindole, ADL5859 exhibited high conformational flexibility and strong interaction with the transmission switch. The putative key residues (W274, D95, V267, L139, V263, M142, T260, R146, R258 and others) involving in the activation pathway were identified through the conventional molecular switch analysis and a pairwise distance analysis, which provides a short list for experimental mutagenesis study. These insights will facilitate further development of therapeutic agents targeting the DOR. Communicated by Ramaswamy H. Sarma


Introduction
Ligands that bind to receptors are typically classified as an agonist or antagonist, meaning they either elicit a biological response or they block the response, respectively. Opioid receptors, members of G-protein-coupled receptors (GPCR), bind an agonist such as opioids outside the cell and adopt the active conformation to activate G protein-dependent signal transduction pathways, ultimately, inducing the biological response of antinociception. These types of proteins are made up of transmembranes (TM) that span the cell membrane (Ritter & Hall, 2009;Shang et al., 2014). There are three types (m, ! and j) of opioid receptors (MOR, DOR and KOR) (Chan et al., 2003), all exhibiting antinociceptive properties. Among them, MOR is the major pharmacological target of current opioid drugs such as morphine. However, these drugs targeting MOR are known for their high propensity for abuse and tolerance, leading to the current opioid epidemic in the United States and the rest of the world. Thus, there is pressing need to develop non-MOR-based opioid drugs (Meert & Vermeirsch, 2005). While drugs acting at the KOR produce dysphoria and anhedonia (Gallantine & Meert, 2008;Knoll & Carlezon, 2010;Land et al., 2008), these adverse effects are not observed when testing drugs that act at the DOR. This makes the DOR a promising pharmacological target for developing new opioid drugs for pain medication (Hudzik et al., 2011).
Over the years, researchers have been able to produce agonists that can target the DOR and cause antinociception and/or antidepressant-like effects without adverse weightbased effects (Gotoh et al., 2017;Nozaki et al., 2014;Saitoh et al., 2013). However, these molecules caused convulsant activities in vivo, making the search for an optimal agonist for the DOR continue. ADL5859, a small molecule agonist targeting the DOR has been discovered with promising properties, including that it does not cause convulsant activities. This ligand completed phase I clinical trials, showing it was well absorbed after oral administration and was generally well tolerated, ready to move on to phase II clinical trials (Corporation, 2008). While this drug has shown promising in vivo activities, there is no crystal structure of ADL5859 with the DOR (Figure 1). Because of this, the detailed interaction of this ligand with the receptor is unknown, as well as its activation mechanism on DOR.
In silico techniques including homology modeling, molecular docking, molecular dynamics (MD) simulations etc. are increasingly popular in studying the interaction between ligands and receptors at atomic level to understand the dynamic behavior of native and mutant receptors (Rajendran, 2016;Rajendran et al., 2012;2016), to understand drug resistance mechanism (Rajendran et al., 2018;Rajendran & Sethumadhavan, 2014), and to identify selective ligands (Kumar Bhardwaj et al., 2021;Sharma et al., 2021). The DOR has not been as well studied as the MOR and KOR. With computational research on the rise, studies were conducted that modeled the DOR even before a crystal structure was available using templates such as bovine rhodopsin. Through this and a 20 ns molecular dynamic simulation, it was suggested that between residues D128 and R192, an internal salt bridge exists. Another study delved into the possible interfaces of DOR homodimeric complexes (Aburi & Smith, 2004;Johnston et al., 2011). After solved crystal structures of the DOR became available, a study that solved two DOR crystal structures themselves modeled the active state DOR and docked small molecules and peptides into their crystal structure to evaluate water-mediated interactions in the binding pocket (Claff et al., 2019). Although the motion of TM6 is a known hallmark of activation for GPCRs, the key change at the residue level is not known for the DOR. Recently, a study (Zhou et al., 2019) has suggested an activation pathway for b 2 -adrenergic receptor, M2 muscarinic receptor, m-opioid receptor, k-opioid receptor, and adenosine A 2A receptor. However, the DOR has not been confirmed to have this pathway, leaving its activation pathway unknown (Zhou et al., 2019).
To give a better insight and understanding of the activation mechanism and detailed interaction of an agonist with the DOR, we ran multiple 1ms molecular dynamic simulations of the novel ADL5859 and of naltrindole starting with a homology model of the inactive DOR (PDB id: 4N6H) by removing the fusion protein and fixing the point mutation (S2P) in the crystal structure. Starting with the inactive conformation and an agonist allows us to probe the conformational change of the receptor toward the active conformation through MD simulations, giving insight into a potential activation mechanism for the DOR. The ADL5859 ligand was first docked to the inactive receptor homology model and this initial structure was used for the simulations. Our most abundant structure from our MD simulation was then compared to the active state crystal complex (PDB id: 6PT3) to find the similarities and differences in the protein-ligand binding interactions and dynamics between the two systems. We observed TM6 move outward with a favorable overlap of our complex with the agonist bound crystal complex, supporting the evidence of the molecular switches for this receptor, specifically the Transmission Switch. Center of mass distances were measured as well using known residue contacts for some class A GPCRs as well as unique distances to our study to determine the conformation change occurring after beginning with the inactive state. The motion of TM6, as well as our additional distance analyses, offers insight into the activation mechanism of the DOR bound with agonist ADL5859. Given the hallmark movements and known residues, we believe our simulation has reached the active-like conformation. Our findings allows for an experimental mutagenesis study using our key residue pathway. The detailed interaction of ADL5859 with the DOR will also help to develop the opioid class of drugs.

Multiple sequence alignment and analysis
Conservation of residues among 62 members of class A GPCR was investigated using Jalview program (Waterhouse et al., 2009). Protein sequences of 61 class A GPCRs listed in the GPCR-EXP database (https://zhanglab.org/GPCR-EXP/) and 5 class A GPCRs listed in the subfamily A4 of rhodopsin family in the Pfam database (http://pfam.xfam.org/family/ PF00001) were retrieved from the UniProt database and MUSCLE tool (Edgar, 2004) with default settings was used for the multiple sequence alignment. Positional occupancy of each residues was analyzed at both family and subfamily A4 levels. The details of sequence conservation are listed in Tables S1 and S2.
Comparison on three crystal structures (4N6H, 6PT2 and 6PT3) of DOR Recently, three crystal structures of the DOR have been solved; one is of the inactive state receptor with an amino terminal b 562 RIL fusion protein in complex with the antagonist, naltrindole (PDB ID: 4N6H), and two structures are of the active state DOR with the peptide agonist KGCHM07 (PDB ID: 6PT2) and with the small molecule agonist DPI-287 (PDB ID: 6PT3), both with a b 562 RIL fusion protein that contained three point mutations. When aligning the crystal structures ( Figure S1), the two agonist bound complexes (6PT2 & 6PT3) show significant overlap, but in contrast, the antagonist complex (4N6H) shows great differences mainly in the fusion protein, as well as TM5, intracellular loop (ICL) 3, and TM6. The difference in TM6 between the crystal structures appears to be one of the molecular switches in GPCR activation known as the Transmission Switch, and is a hallmark of GPCR activation, indicating the agonist bound complexes as the active crystal structures and the antagonist bound receptor as the inactive crystal structure. However, because the fusion protein also adopts very different conformations in the two complexes, the probability of the fusion protein causing the conformational change of TM6 cannot be completely ruled out. Because both DPI-287 and ADL5859 share the same pharmacophore and biological activity, their action on the DOR should be similar ( Figure 1).

The homology modeling of the human DOR in an inactive conformation
The FASTA files of human DOR-1(P41143) were obtained from the UniProt website (Magrane, 2011) to fix any nonnative mutations or fusion protein in the homology modeling using the inactive DOR conformation (PDB ID: 4N6H) as the template structure. This inactive DOR conformation allows us to investigate the activation mechanism of DOR by agonist ADL5859. The sequence comparison between the native sequence and the sequence in 4N6H is shown in Figure S5 (Granier et al., 2012;Huang et al., 2015). The Prime feature in Maestro (Jacobson et al., 2002; was used to build the homology model based on pre-aligned crystal structure of human d-opioid 7TM receptor bound with antagonist naltrindole (PDB ID: 4N6H) oriented in the membrane from the OPM web server site (Lomize et al., 2012). The crystal structure of human d-opioid 7TM receptor bound with agonist DPI-287 (PDB ID: 6PT3) was used for comparison. Alignment is based on the backbone residues of TM1-5&7 (4N6H: .

Molecular docking of ADL5859 into the inactive DOR conformation
Protein structures were prepared using Maestro's Protein Preparation Wizard (Sastry et al., 2013). The charge state of the preprocessed protein was optimized at pH 7. A restrained minimization was then performed to relax the protein structure using OPLS3 force field (Harder et al., 2016). The 3D structures of naltrindole and ADL5859 were prepared using Maestro Elements. Figure 1 shows the chemical structures of the ligands. The 2D structure of ADL5859 was drawn using Maestro's 2D sketcher, converting into 3D structure. 3D structure of naltrindole was extracted from the crystal structure (PDB ID: 4N6H). Both ligands had their bond orders fixed and underwent the pKa calculation. The ionization/ tautomeric states were generated at pH 7 using Maestro's Epik tool based on the more accurate Hammet and Taft methodologies (Sastry et al., 2013). The lowest ionization/ tautomeric state (þ1) was selected and the geometry was then minimized to the most energetically favorable structure to relax the structure of the ligand. Because the inactive DOR has been crystallized, the binding site is known. The receptor docking grid files were generated from the prepared receptor, in which the centroid of the crystal ligand, naltrindole, was used to specify the active site. Since the crystal ligand (naltrindole) is bound to the orthosteric site, the agonist ADL5859 should also bind to this pocket. The prepared ligands (naltrindole and ADL5859) were docked into their corresponding generated grids using Glide XP scoring with default procedures and parameters . In detail, the receptor grid required for the docking process was generated using van der Waals scaling factor of 1 and partial charge cutoff 0.25. Docking was performed using a ligand-centered grid using OPLS3 force field (Harder et al., 2016). Glide XP Dock performed a comprehensive systematic search for the best receptor conformations and orientations to fit the ligand. The docked poses were compared to the inactive crystal complex (PDB ID: 4N6H) with an antagonist to verify if the docked ligand poses were reasonable. Both ligands were bound within the binding pocket with naltrindole binding similar to the crystal ligand with subtle differences, providing a reasonable starting pose for later molecular dynamic simulations. The binding pose can then be refined given the full conformation flexibility in the simulations.

Molecular dynamic simulation system setup
Two MD simulation systems were constructed using the prepared and refined receptor-ligand complex from the Glide XP docking and crystal structure (4N6H) as input files. Each system was built using SPC as water solvent model (Mark & Nilsson, 2001) using orthorhombic solvent box with 10 Å water buffer. System was neutralized using Na þ and Clions and were added to maintain the salt concentration of 0.15 M NaCl. After the system was successfully solvated, OPLS3 force field (Harder et al., 2016) was used to represent the receptor-ligand.
Using Desmond module, the system was first relaxed using the default relaxation protocol for membrane proteins (Zhang et al., 2012). This relaxation protocol consists of eight stages: 1). Minimization with restraints on solute heavy atoms; 2). Minimization without any restraints; 3). Simulation with heating from 0 K to 300 K, H 2 O barrier and gradual restraining; 4). Simulation under the NPT ensemble (constant number of particles, constant pressure of 1 bar and constant temperature of 300 K) with H 2 O barrier and with heavy atoms restrained; 5). Simulation under the NPT ensemble with equilibration of solvent and lipids; 6). Simulation under the NPT ensemble with protein heavy atoms annealing from 10.0 kcal/mol to 2.0 kcal/mol; 7). Simulation under the NPT ensemble with Ca atoms restrained at 2 kcal/mol; and 8). Simulation for 1.5 ns under the NPT ensemble with no restraints. After the relaxation, three trajectories (1000.0 ns of each) were conducted under the NPT ensemble for each of the systems using the default protocol. In details, temperature was controlled by using the Nos e-Hoover chain coupling scheme (Ikeguchi, 2004) with a coupling constant of 1.0 ps. Pressure was controlled using the Martyna-Tuckerman-Klein chain coupling scheme (Ikeguchi, 2004) with coupling constant of 2.0 ps. M-SHAKE (Bailey & Lowe, 2009) was applied to constrain all bonds connecting hydrogen atoms, enabling a 2.0 fs time step in the simulation. The k-space Gaussian split Ewald method (Shan et al., 2005) was used to treat long-range electrostatic interactions under periodic boundary conditions (charge grid spacing of $1.0 Å, and direct sum tolerance of 10 À9 ). The cutoff distance for short-range nonbonded interactions was 9 Å, with the long-range Van der Waals interactions based on a uniform density approximation. To reduce the computation, non-bonded forces were calculated using an r-RESPA integrator (Stuart et al., 1996) where the short-range forces were updated every step and the long-range forces updated every three steps. The trajectories were saved at 40.0 ps intervals for analysis.
To check the convergence of MD simulations, the protein Ca and ligand RMSD plots were investigated for each trajectory. The relatively flat plots within last 200 ns indicates that the complex systems have reached a steady state.
Desmond SID tool was used to analyze the behavior and interaction of proteins and ligands during the course of simulation including RMSD; protein-ligand contacts including H-bonds, hydrophobic, ionic, and water-bridge contacts; secondary structure changes and Root Mean Square Fluctuation (RMSF) measures.

Trajectory clustering analysis
Desmond trajectory clustering tool (Bowers et al., 2006) was used to group complex structures from the last 100 ns simulation of each complex system. Backbone RMSD matrix was used as structural similarity metric, the hierarchical clustering with average linkage (Bowers et al., 2006) was selected as the clustering method. The merging distance cutoff was set to be 2 Å. The centroid structure (i.e. the structure having the largest number of neighbors in the structural family) was used to represent the structural family. The centroid structures of populated structural families (>2% of total structure population) are shown in Figures S10-S12 of the supporting material.

RMSD of protein structure analysis
We analyzed the RMSD of the whole protein, transmembrane domain (TMD), and TM6 for each trajectory in reference to both the solved inactive crystal structure of the DOR (PDB ID: 4N6H) and the solved active crystal structure of the DOR (PDB ID: 6PT3) to aid in differentiating the protein in the active or inactive state. Alignment of the whole protein is based on the backbone residues of the protein excluding backbone residues of TM6 (4N6H: 36-248, 288-321, 6PT3: 45-248, 288-321), and alignment of the 7TMs is based on backbone residues of TM1-5&7 (4N6H: , both of which exclude oxygen. RMSD calculations are of the backbone residues of the whole protein (4N6H: 36-321, 6PT3: 45-321), 7TMs (4N6H: , and TM6 (250-286). Maestro's Simulation Event Analysis was used to calculate the RMSD of the MD simulation. VMD RMSD Calculator (Humphrey et al., 1996) was used to calculate the RMSD of the portions of each representative structure of abundant clusters from each trajectory.

Analysis of microswitches and pairwise residue distance
We analyzed three important microswitches (Trzaskowski et al., 2012), i.e. transmission (CWXP), tyrosine toggle (NPXXY) and ionic lock (DRY) which are used to differentiate between active and inactive states of class A GPCR. Another lock known as 3-7 lock (Trzaskowski et al., 2012) was also investigated.
Framewise center of mass distances were calculated for 23 pairs of residues. Out of 23 pairs of residues, 20 were previously used by Zhou et al. (Zhou et al., 2019), to define the activation of mechanism of class A GPCRs. For calculation of atomic and center of mass distances, VMD (Humphrey et al., 1996) was used.

Normal mode analysis
The individual trajectories from both systems were used in the Normal Mode Wizard in VMD (Bakan et al., 2011) to generate a principal component analysis of the top 10 normal modes.

Sequence-based comparison of human DOR with other class A GPCRs
Comparison of primary sequence of human DOR with 61 class A GPCRs revealed that DOR sequence was more conserved at subfamily level than the family level (Tables S1-S4 and Figures S3 and S4). Residues belonging to molecular switches such as CPXY, NPXXY and D/ERY motifs showed high level of conservation at both family and subfamily levels. DOR showed more than 50% identity with kappa (j) and mu (m) opioid receptors.

Homology model of human DOR without fusion protein
Homology models of the DOR for the generation of ADL5859 and naltrindole bound complexes were obtained based on the experimentally resolved inactive crystal structure of the human DOR with naltrindole (PDB ID: 4N6H) as described in the methods section. The inactive crystal structure of the DOR was used as a template with the full length human DOR sequence filling in missing side-chain information and fixing any non-native mutations. In the homology model of the DOR, mutation P37S present in the crystal structure of DOR was corrected. A cartoon representation is shown in Figure S5 from top and side viewpoints. This homology model is a valid starting structure for molecular docking and simulations to inspect the activation mechanism.
Docked conformation of ADL5859 shows similarity with the co-crystallized conformation of DPI-287 and subtle difference with naltrindole The docking protocol was validated by docking the co-crystallized naltrindole back into the crystal structure of DOR with its crystal pose ( Figure 2C). The result had an almost identical overlap which suggests a satisfactory docking method. ADL5859 was docked into the inactive receptor, resulting in a similar binding pose as the crystal agonist ( Figure 2B). ADL58589 and naltrindole do not show overlap as they have different scaffolds but maintain the same space in the binding pocket ( Figure 2A). Docked conformation of ADL5859 showed interaction with salt bridge with D128 and non-polar interactions with M132, W274, V281, L300 and Y308.

MD Simulations of human DOR with ADL5859 and naltrindole assume steady state
The complexes from Glide XP docking were used as the input structures to perform six independent 1 ms molecular dynamic simulation for ADL5859-DOR complex and three 1 ms molecular dynamic simulation for naltrindole-DOR complex. The use of multiple simulations adds to the statistical significance of the simulation results. The RMSD values of ADL5859 and naltrindole shows the flexibility of the ligand interaction with the DOR. The movement of the ligand and the receptor during the simulation was analyzed using the RMSD (Figure 3 and S6 and S7). In the RMSD plots for systems, the protein and ligand RMSD are approaching stable values within the last 200 ns, indicating that the two complexes had been sufficiently equilibrated and allowed for rigorous analysis. The protein RMSD for ADL5859 changes $3Å which indicates that the protein is undergoing a large conformational change during the simulation ( Figure 3). The protein RMSD for naltrindole has a lower RMSD value and remains more stable. The ligand RMSD for ADL5859 shows multiple fluctuations around 200-400 ns during the simulation ( Figure 3). Naltrindole RMSD has a large fluctuation around 200 ns and again between 800 and 950 ns. ADL5859-DOR system shows adopting toward the crystal active conformation while naltrindole-DOR system maintains inactive conformation To determine whether our agonist MD complex can adopt toward the active conformation and move away from the inactive conformation, the RMSD of the receptor against the active and inactive crystal structures was calculated throughout the whole trajectory. The RMSD of the whole protein, TMD, and TM6 of the receptor was calculated to home in on which portion of the protein has the most conformational change, presumably leading to activation. The RMSD of the same regions for the antagonist naltrindole system was calculated as well to determine if the MD complex will maintain the inactive conformation to serve as a control. For this to be clear, our complexes were aligned against the inactive and active crystal structures (4N6H & 6PT3) using the backbone residues of the TMD excluding TM6. The reason for this exclusion is because TM6 appears to have the most dynamics while the other transmembranes are more rigid. Aligning . Protein (blue) and ligand (red) RMSD averaged over the six independent trajectories for agonist ADL5859-DOR system (A) and three independent trajectories for antagonist Naltrindole-DOR system (B). Alignment is based on the Ca atoms of the protein and the RMSD calculation with reference to the first frame. Ligand RMSD calculated the ligand heavy atoms after the complex is first aligned on the protein backbone of the reference. The data for individual trajectories can be found in Figures S6 and S7 of the supporting document.
the rigid parts of the receptor will better show the change in conformation of the flexible part. Examining the whole protein of the individual agonist RMSD plots ( Figure S8) revealed that trajectories 1, 3, 5, and 6 behave similarly as they clearly show a general trend where they increase in RMSD when against the inactive crystal structure and decrease in RMSD when against the active crystal structure. This trend becomes more pronounced when looking at the RMSD of the TMD for these trajectories and even more so for the RMSD of TM6. Using trajectory 1 as an example, the RMSD value of the TMD at the beginning of the trajectory against the inactive crystal structure is 1.9 Å and ends at 2.4 Å; and when against the active crystal structure, the RMSD value begins at 3.0 Å and ends at 2.4 Å. With this same process for TM6 in trajectory 1, the RMSD values when it first begins against the inactive crystal structure is 3.3 Å and ends with an RMSD value of 4.2 Å; and when against the active crystal structure, it begins with an RMSD value of 5.7 Å and ends with an RMSD value of 4.0 Å. This indicates the major part contributing to the protein folding is the movement in TM6. Trajectory 4 of the agonist system shows the same trend within the first 800 ns as trajectories 1, 3, 5, and 6 but after 800 ns, it reverses the process with an increasing trend against the active crystal structure and decreasing trend against the inactive crystal structure. Trajectory 2 of the agonist system differs from the others in which it only maintains the same RMSD trend as the other trajectories during the first 200 ns but still not as clear. After 200 ns, especially for the RMSD of TM6, it loses the general trend and levels out. To clarify what this means, using the RMSD of TM6 in trajectory 2, when against the inactive conformation at 200 ns has an RMSD value of 2.3 Å and ends with an RMSD value of 2.4 Å; and when it is against the active conformation at 200 ns has an RMSD value of 4.3 Å and ends with an RMSD value of 4.5 Å. In the case of the RMSD for the antagonist naltrindole system ( Figure S9), all three trajectories displayed the same pattern for each protein region in which there was no change in RMSD when against either crystal structures. From these results we can speculate that if the system adopts toward the active conformation, it should decrease in RMSD against the active crystal, favoring smaller RMSD values, and increase in RMSD against the inactive crystal, favoring larger RMSD values. If there is no protein folding within the system, then the RMSD should not change which would be validated by flat RMSD plots. Because the majority of the agonist trajectories shared a general trend in RMSD and all antagonist trajectories shared a general trend in RMSD, the individual trajectories were averaged to simplify the results (Figure 4). After averaging the trajectories for both systems, a noticeable feature in the RMSD became more apparent. For both simulation systems (agonist & antagonist), there is a gap between the starting values of the RMSD when against the reference structures, with the inactive reference structure at much lower values than when against the active crystal structure. Our simulation systems began using the inactive conformation which correlates with this feature that our systems are more similar to the inactive crystal structure at the beginning of the trajectory, depicted by lower RMSD values. This gap gets larger as the protein regions get more specific (whole protein ! TMD ! TM6), indicating the difference between the conformation of the reference structures increases as the RMSD gets more specific. The averaged agonist trajectories show to have maintained the same clear trend in RMSD as its individual trajectories throughout the protein regions, with the gap diminishing between the reference structures throughout the Figure 5. MD DOR-agonist-ADL5859 structure (purple&blue) of the most abundant cluster (52%) of the first trajectory, its comparison to the solved crystal structures (Inactive PDB ID: 4N6H, cyan&red) (Active PDB ID: 6PT3, tan&yellow) in the side, top, and bottom views respectively. RMSD of the protein, 7TMs, and TM6 of the MD structure against crystal structures is listed. Alignment is using backbone residues of TM1-5&7, excluding oxygen. Ligand view of MD-DOR-ADL5859, crystal antagonist Naltrindole in inactive DOR, and crystal agonist DPI-287 in active DOR. The other clusters from trj1 and the clusters from the other trajectories can be found in Figure S10 of the supporting document.
trajectory, and the most prominent trend occurring in TM6. RMSD values of the different protein regions at the beginning and end of the simulation for the averaged simulation systems are clearly tabulated in Table S3 of the supporting document. The RMSD of the averaged agonist system of the whole protein against the active crystal begins with an RMSD of 3.5 Å and ends with a decreasing RMSD of 3.4 Å showing a difference of .1 Å. When against the inactive crystal, the whole protein begins with an RMSD of 2.4 Å and ends with an increasing RMSD of 3.4 Å, showing a difference  of 1 Å. The TMD displays the trend in RMSD more clearly with the average agonist trajectories against the active crystal beginning at an RMSD of 2.9 Å and ending with an RMSD of 2.7 Å, showing a decreasing difference of 0.2 Å. When against the inactive crystal, the TMD begins with an RMSD of 1.6 Å and ends with an RMSD of 2.6 Å, showing an increasing difference of 1 Å. The greatest difference in RMSD lies in TM6 beginning with an RMSD of 5.6 Å against the active crystal and ending with a decreasing RMSD of 4.4 Å, indicating a 1.2 Å difference. When against the inactive crystal structure, TM6 of the averaged agonist trajectories begins with an RMSD of 2.6 Å and ends with an increasing RMSD of 3.6 Å, showing a 1 Å difference. The trend in RMSD of the whole protein and the TMD of the agonist system was not as drastic, indicating TM6 shows the most change in conformation. When viewing the RMSD of the averaged antagonist trajectories, they have still maintained the trend of no change in RMSD with relatively flat plots. The gap between the reference structures is maintained throughout the whole trajectory for each protein region, indicating the antagonist system is more similar to the inactive crystal structure and less similar to the conformation of the active crystal structure. When against the active crystal structure, the averaged antagonist system of the whole protein begins with an RMSD of 3.1 Å and ends with an RMSD of 3.2 Å, showing a .1 Å difference. When against the inactive crystal structure, the RMSD of the whole protein begins at 3.3 Å and ends with an RMSD of 3.3 Å, indicating no change from the inactive reference structure. The RMSD of the TMD when against the active crystal structure begins at 2.9 Å and ends with an RMSD of 3.0 Å, showing an increase of .1 Å. When against the inactive crystal structure, the TMD starts with an RMSD of 2.4 Å and ends with an RMSD of 2.4 Å, again showing the antagonist system has not folded away from the inactive conformation. The RMSD of TM6 of the averaged antagonist system when against the active crystal structure begins at 5.2 Å and ends at an RMSD of 5.3 Å, showing a 0.1 Å increase. When against the inactive crystal structure, TM6 begins with an RMSD of 2.4 Å and ends with an RMSD of 2.5 Å, showing another 0.1 Å increase. This data alludes to our agonist system folding towards the active conformation as the gap between the reference structures dies out with a decreasing trend against the active crystal structure and an increasing trend against the inactive crystal structure. Lower RMSD values against the active crystal structure indicate greater similarity with the active conformation and in this same respect, higher values against the inactive crystal structure indicate the difference from the inactive conformation. This trend of lower values against the active DOR structure (6PT3) and higher values against the inactive DOR structure (4N6H) was the most prominent for the RMSD of TM6 of the agonist system, in which the movement of TM6 is associated with activation in GPCRs, specifically the Transmission Switch as previously mentioned in the introduction. In the case of the antagonist system, the RMSD of the different protein regions showed either a change of 0.1 Å or no change at all while also maintaining the gap between reference structures, which would imply that the system has not folded away from the inactive conformation. From the backbone RMSD plots, it is clear the activation related conformational changes are unique to the agonist-ADL5859 bound DOR.

ADL5859-DOR system has multiple clusters folding toward the active state
The clustering analysis was done after the MD simulation to identify the populated conformation for each trajectory. Based on the clustering algorithm, each cluster conformation contains a percentage of abundance in which a cutoff of 1% was used. The clustering results are clearly tabulated in Table S5. There was a total of 15 abundant clusters from the six agonist trajectories with each trajectory containing at least two abundant conformations, with the exception of trajectory 2. For the three antagonist trajectories, there were three abundant clusters, one for each independent trajectory. The clustering result for both systems corresponds with the thought that the agonist system is believed to be folding towards the active state which would produce more than one conformation; and in this same regard, the antagonist system should not have any folding, which would result in one conformation.
To acquire a deeper insight into the conformation of the MD structures, the RMSD of the whole protein, TMD, and TM6 of each abundant cluster against the inactive and active crystal structures was calculated as described in the Methods  (Table S5 and Figures S10-S12). To consider an abundant cluster as folded toward the active state, when comparing the RMSD values of the protein regions against the different crystal structures, especially in the case of TM6, they should differ roughly by at least 1 Å, with lower values when against the active crystal structure, indicating greater similarity. The  2.0 ± 0.9 1.6 ± 0.7 2.2 ± 0.9 0.6 0.9 ± 0.3 0.8 ± 0.0 1.1 ± 0.2 0.2 3 Y318 7.53 -Y318 7.53 1.9 ± 0.5 1.7 ± 0.4 1.9 ± 0.3 0.2 1.0 ± 0.4 0.9 ± 0.1 1.0 ± 0.4 0.1 4 R146 3.50 -T260 6.34 5.2 ± 1.6 4.6 ± 0.1 6.5 ± 2.6 1.9 4.6 ± 0.4 4.6 ± 0.0 4.6 ± 0.1 0.0 5 I226 5.51 -F270 6.44 11.4 ± 0.7 11.7 ± 0.4 11.1 ± 0.6 À0.6 11.4 ± 0.5 11.5 ± 0.0 11.2 ± 0.1 À0.3 6 I136 3.40 -W274 6.48 7.5 ± 0.6 7.4 ± 0.7 7.7 ± 0.3 0.3 7.4 ± 0.4 7.3 ± 0.2 7.6 ± 0.1 0.3 7 T230 5.55 -V267 6.41 10.4 ± 0.6 10.1 ± 0.4 10.5 ± 0.2 0.4 9.5 ± 0.6 9.6 ± 0.0 9.4 ± 0.5 À0.2 8 L139 3.47 -N314 7.49 11.2 ± 0.5 11.1 ± 0.4 11.1 ± 0.3 0.0 11.6 ± 0.6 11.8 ± 0.4 11.3 ± 0.3 À0.5 9 Y233 5.58 -V266 6.40 12.0 ± 1.2 11.7 ± 1.2 11.6 ± 0.9 À0.1 10.1 ± 0.6 9.8 ± 0.0 10.5 ± 0.4 0.5 10 L139 3.43 -Y318 7.53 11.5 ± 0.6 11.7 ± 0.2 11.3 ± 0.3 À0.4 11.9 ± 0.5 12.4 ± 0.0 11.7 ± 0. RMSD values of the abundant clusters from the antagonist system against both crystal structures were also taken into consideration when determining the criteria for folding toward the active state as this system should remain in the inactive state. Each cluster of the antagonist system had an RMSD ranging from 3.1 to 3.4 Å for the whole protein when against the active crystal structure and 3.1-3.3 Å when against the inactive crystal structure. When looking at the TMD, the antagonist abundant clusters had an RMSD ranging from 2.8 to 3.0 Å when against the active crystal structure and 2.7-2.8 Å when against the inactive crystal structure. The last protein region calculated for the RMSD of the naltrindole complexes, TM6, had RMSD values ranging from 4.9 to 5.5 Å against the active crystal structure and 1.8-2.9 Å against the inactive crystal structure. After calculating the RMSD of each of the agonist clusters against the crystal structures, there were three trajectories (trj1, trj4, trj6) and four clusters (trj1 cluster 1&3, trj4 cluster 5, trj6 cluster 1) that can be considered in the active state. The most abundant cluster (cluster 1) from trajectory 1 in the agonist system is the most favorable from this system and is used as the representative structure ( Figure 5). The interactions of agonist ADL5859 with the DOR from the first trajectory were examined, along with the crystal antagonist naltrindole in the crystal inactive state DOR, and the crystal agonist DPI-287 in the crystal active state DOR. With this, ADL5859 directly interacted with polar histidine 278 of the DOR during the simulation, which is consistent with the interactions of crystal agonist DPI-287 and the DOR. Both agonists also have exposure to solvent on the same CH 3 of their structures. The clusters from the naltrindole system show good agreement with the inactive crystal structure and were validated by the RMSD values with the abundant cluster from the first trajectory used as the representative structure (Table S5, Figures S11 and S12). Naltrindole directly interacted with negatively charged aspartic acid 128 of the DOR during the simulation which is consistent with the interactions of crystal antagonist naltrindole and the DOR. This indicates the conformational changes occurring to the receptor are unique to the agonist ADL5859-DOR system. While the RMSD of the superimposed agonist structures indicated four clusters folding toward the active state, there were other abundant clusters from the agonist-ADL5859 system that are thought to contain intermediate conformations of the protein folding toward the active state ( Figure S10). These other noteworthy structures are cluster 1 and cluster 2 from trajectory 3 with a percentage of abundance of 90 and 11 respectively; trajectory 4 cluster 1, cluster 3, and cluster 4 with percentages of abundance of 43, 20, and 5 respectively; and cluster 1 from trajectory 5 with a percentage of abundance of 56. Even though these structures did not have as favorable of RMSD values as the others that were identified, they are speculated to be intermediate conformations in the activation-related folding process, as they also still display agreement with the active crystal structure and contain great differences from the inactive crystal structure. The differences from the inactive crystal structure mentioned mainly lie in TM5, ICL3, and TM6 while similarities to the active crystal structure mainly lie in TM1, TM5, and TM6. These clusters, along with the ones identified from the RMSD as folding toward the active state, are consistent with the results from the backbone RMSD plots ( Figures S8 and S9), in which the trajectories these cluster conformations were generated from (trj1, trj3, trj4, trj5, trj6) displayed the general trend implying protein folding toward the active conformation throughout the simulation. This adds further support of the agonist-ADL5859 system folding toward the active state.
Secondary structure analysis shows maintained helices but with subtle differences between the two complexes Secondary structure elements (SSE) plots illustrate where the major differences occur in the transmembrane sections for comparisons to be made between the agonist and antagonist ( Figure S13). Some of the noteworthy features include: A) Additional kink in ADL5859 in TM6. B) Additional kinks in naltrindole in TM2 and TM7. SSE plots for the individual ADL5859 and naltrindole trajectories can be found in Figures S14 and S15, respectively.

Key interactions of ADL5859 and naltrindole with human DOR
The Desmond Simulation Interaction Diagram was used, as described in the methods section, to analyze and compare the residues involved in both naltrindole and ADL5859 binding to the receptor ( Figure 6). The key interacting residues that interacted with each ligand of the combined trajectories for more than 5% of the simulation are tabulated for clarity in Table S4. Key interacting residues of the individual trajectories for the agonist and antagonist systems that lasted more than 20% of the MD simulation can be found in Figures S16 and S17, respectively. Interacting residues during the simulation in histogram format of the individual agonist and antagonist trajectories can be found in Figures S18 and S19, respectively. Using both Figure 6 and Table S4 it can be understood that ADL5859 interacted with a greater number of residues from the receptor with a total of 47 interacting residues that mainly occur in TM2, TM3, TM5, TM6, and TM7. Naltrindole had a total of 30 interacting residues in which 26 of them were conserved in ADL5859. In TM6, there were 6 conserved interacting residues between the two systems. Naltrindole, however, interacted with residue I282 6.56 which differed from ADL5859 that interacted with residues F270 6.44 and W285 6.59 . The major differences lie in TM7 where ADL5859 had 10 interacting residues: V297 7.32 ; L300 7.35 ; H301 7.36 ; C303 7.38 ; I304 7.39 ; G307 7.42 ; Y308 7.43 ; N310 7.45 ; S311 7.46 ; N314 7.49 whereas naltrindole had only 6: L300 7.35 ; I304 7.39 ; G307 7.42 ; Y308 7.43 ; N310 7.45 ; and S311 7.46 . The 2D interaction diagram indicates that p-p stacking is crucial for the binding of both ligands ( Figure 6).

DOR shows higher flexibility when bound to ADL5859 than naltrindole
The receptor RMSF (Figures 7 and S20-S21) shows higher RMSF values in areas of high flexibility such as the intra-and extracellular loops, as well as the N-and C-terminals. In contrast, the helices, or rigid areas of the receptor, show lower RMSF values (Table 1). When the DOR was bound to ADL5859, the receptor RMSF values were higher than when bound to naltrindole. ADL5859 has higher ligand RMSF values in comparison to naltrindole (Figure 7). In Figure 7, naltrindole has lower RMSF values at almost every atom. This is consistent with what was expected from the structures of the two ligandsfor naltrindole contains one more six membered ring than ADL5859, as well as two five membered rings and one three membered ring, making naltrindole a more rigid ligand compared to ADL5859. These differences indicate that the DOR exhibited greater flexibility when interacting with an agonist (ADL5859) than with an antagonist (naltrindole). This is likely because the receptor was undergoing activation when in complex with ADL5859. Receptor RMSF data of the individual trajectories can be found in Figures S20 and S21 of the supporting document and individual ligand RMSF data can be found in Figures S22 and  S23 of the supporting document.
Distance of ADL5859 with side chain of D95 decreases while naltrindole shows no change Residue D95 2.50 plays an important role in DOR activation signaling. Mutagenesis studies have shown that G protein signaling is affected in GPCRs when D 2.50 is replaced with a neutral amino acid (White et al., 2018). Measuring the distance of the positively charged nitrogen atom of agonist ADL5859 and antagonist naltrindole with D95 of the protein shows a clear difference between the two systems ( Figure 8). In the agonist ADL5859 system, the distance decreases about 4 Å inferring a salt bridge formed with the ligand. The antagonist system shows no change in distance, further supporting the DOR in complex with the antagonist naltrindole maintained the inactive conformation.

Role of molecular switches in the activation of DOR by ADL5859
Class A GPCRs share a set of conserved motifs which are known as molecular switches. Activation of the class A GPCR can be characterized by the conformational changes in one or more molecular switches. Molecular switches were evaluated throughout the entire trajectories of ADL5859-DOR and naltrindole-DOR systems. For the DOR, three such motifs are considered to be critical in regulating its activity: Transmission Switch (CWXP) in TM6, Ionic Lock Switch (DRY) in TM3 and TM6, and Tyrosine-Toggle Switch (NPXXY, X represents any amino acid) in TM7. Each abundant cluster from the individual trajectories of both simulation systems was compared to the solved crystal structure of the inactive DOR in complex with antagonist naltrindole (PDB ID: 4N6H) to measure the switch distances. We have determined five abundant clusters in the agonist system (trj1 cluster 1&3, trj3 cluster 2, trj4 cluster 5, trj6 cluster 1) to have broken molecular switches, inferring activation, based on the transmission switch and ionic lock switch as there was no significant change in the tyrosine toggle switch for any of the clusters (Figure 9, Figures S24-S31). Trajectory 1 cluster 1 is the most favorable from the agonist system as is used as the representative structure (Figure 9). The key distances characterizing the switches for all trajectories of ADL5859 and naltrindole are summarized in Table S6. The transmission switch is the most observable switch out of the three that when in the active state, a large outward 'swing' occurs in TM6. Our molecular switch analysis shows this same 'swing' in our agonist ADL5859 system ranging from 9.2 Å to 13.9 Å. The naltrindole system displays a motion of 3.8 Å and 3.0 Å in TM6 indicating it lacks any significant outward motion of TM6. One of the aromatic rings of ADL5859 showed interaction with the side chain of W274, a highly conserved residue in the transmission switch. This is consistent with the active crystal structures of the DOR exhibiting slight rotameric change at the level of W274. Naltrindole also interacts with W274 but interaction is weak compared to ADL5859. The salt bridge between D 3.49 and R 3.50 residues in the DRY motif is already disrupted in the inactive DOR with distances of >3.5 Å and is missing an acidic residue in position 6.30, which is only conserved in about 30% of GPCRs. To compensate for this, we measured the hydrogen bonding between residues R146 3.50 and T260 6.34 , for this interaction is thought to stabilize receptors in the inactive state and may be important for regulating receptor signaling (Katritch et al., 2013). When measuring the distance of the agonist-ADL5859 system, the hydrogen bonding between residues R146 3.50 and T260 6.34 was disrupted in trj1 cluster 1, trj3 cluster 2, trj4 cluster 5, and trj6 cluster 1 with distances ranging from 4.2 Å to 9.3 Å, implying the inactive state of the receptor was destabilized. The distances for the ionic lock switch in the antagonist-naltrindole system ranged from 3.7 to 4.4 Å, indicating it was only broken in one trajectory (trj 2).

Center of mass distances identify activation related changes
The conserved motifs in molecular switches of Class A GPCRs were used as the motivation to search for other potential motifs in examining the conformational changes throughout the simulation systems (Table 2). From the literature, 18 pairwise residues were identified that were analyzed in this study, including several hydrophobic and van der Waals contacts, as well as hydrogen bonding and salt bridges (Zhou et al., 2019). In addition, four atomic mass distances were measured from the molecular switch motifs: S249-S249 in TM6, W274-W274 in TM6, Y318-Y318 in TM7, R146-T260 hydrogen bonding between TM3-TM6 (Claff et al., 2019). When a residue is a part of TM6, the distance between that and the other transmembrane helices, particularly TM3 and TM7, should increase, breaking any interactions. Based on this, four pair-wise residues involving TM6 that are also conserved in the subfamily (opioid receptors) level were additionally measured: V265-L321 Van der Waals contact, R258-E323 salt bridge, L139-V266 Van der Waals contact, I88-M262 hydrophobic and Van der Waals contact. In total, we analyzed 27 distances, measuring the center of mass between residues, to investigate the conformational changes over the whole trajectory for the average of both simulation systems (Table 3). Through this, seven distance profiles were identified that are associated with activation. A significant increase in the distance of the center of mass is observed in the agonist system when one of the residues in a pair is a part of TM6. The timelines for the distances involving residues in TM6 for both systems are shown in Figure S32. In the agonist-ADL5859 system, trajectory 1 and 6 displayed similar distance profiles with a general increasing trend, and because of this, were averaged together to simplify the distance timeline. The three antagonist-naltrindole trajectories displayed similar distance profiles as well with a general slightly decreasing trend, and therefore, were averaged together to simplify the distance timeline. The greatest difference between the two simulation systems can be observed between residues I88 2.43 and M262 6.36 located at the intracellular portions of TM2 and TM6. Over the course of the simulation for the agonist system, this residue shows a large increase throughout the ADL5859 timeline with a change of 4.1 Å. The timeline of the naltrindole system shows this residue to be gradually decreasing throughout the simulation with a change of 0.2 Å. This is consistent with the previously mentioned thought that the distance between TM6 and other transmembrane helices should increase.

Normal mode analysis shows the overall motion of the receptor
The normal mode analysis ( Figures S33 and S34) shows the top 10 low vibrational modes identified from the principal component analysis, however modes 1-3 show the most difference in motion. The overall motion of the receptor is significant, however, the motion of TM6 is presumed to be a hallmark in determining the state of the receptor (i.e. active, inactive) and therefore can be focused in on for clarity. Using our representative structures from our two systems, our analysis identified the ADL5859 system having a higher degree of fluctuation and a more distinct motion when compared to naltrindole. In the top mode of the ADL5859 system, the ends of TM5, TM6, and ICL3 appear to be moving outward, deeming the bottom of the receptor to be open, or otherwise, in the active state. The motion in the naltrindole system of the same region can be observed as moving back towards TM3 and TM4, still in a closed position. This analysis was also performed on trajectory six of the ADL5859 system and the remaining two trajectories of the naltrindole system. Trajectory six of the ADL5859 system had maintained similar modes of motion from the first agonist trajectory-. The two trajectories from the naltrindole system also maintained similar modes of motion from the first antagonist trajectory.

Discussion
Computational studies have been utilized in studying the DOR even before crystal structures became available, however, to our best knowledge, no known studies have applied the use of MD simulations to study its activation mechanism started from an inactive conformation. Molecular dynamic simulations is extremely powerful tool to probe deeper into the interactions and dynamics occurring in a protein system. Starting from the inactive state of the DOR and sampling its active conformation in the presence of a small molecule agonist offers an opportunity to explore the structural basis of DOR activation. Using agonist ADL5859 adds to the novelty due to its lack of a crystal structure and unknown interactions with the DOR. ADL5859 shares a pharmacophore with the crystal agonist DPI-287 (Figure 1), implying their interactions with the DOR will be similar, validating the use of this small molecule ligand for this study. The crystal antagonist naltrindole was used as our second simulation system in order to determine if the inactive conformation of the DOR when bound to an antagonist can be maintained in MD simulations to serve as a negative control. Here we present the first study probing the active conformation of the DOR started from an inactive conformation using MD simulations.
After docking the ligands into the homology model, both agonist ADL5859 and antagonist naltrindole resulted in docking conformations consistent with the crystal ligand poses ( Figure 2) (Claff et al., 2019;Fenalti et al., 2014). This validated our initial poses for the MD simulation of both ligands. Following the simulation of both systems, the calculation of RMSD throughout the simulation of different protein regions against both crystal structures (4N6H, 6PT3) offered a deeper insight into the structural change response of the DOR when bound to an agonist or antagonist (Table S3, Figure 4). The gap between the crystal reference structures in the RMSD plots indicated both simulation systems were more similar to the inactive conformation at the start of the simulation, which is consistent with our starting structures that were generated based on the inactive crystal structure (4N6H). Decreasing in RMSD against the active crystal structure while increasing in RMSD against the inactive crystal structure would suggest that the system is adopting toward the active conformation. Encouragingly, the agonist-ADL5859 system followed this RMSD trend with the most apparent in TM6. This correlates with what is known about GPCRs, that TM6 outward opening is a hallmark of the activation in class A Figure 10. Binding of agonist to DOR led to flow of allosteric signal from layer 1 (containing agonist binding residues) to different layers of human DOR.
GPCRs (Trzaskowski et al., 2012). The RMSDs in the antagonist system only showed no change. This analysis supports our agonist system folding toward the active conformation and our antagonist system maintaining the inactive conformation.
The trajectory clustering analysis identified the populated conformations of the individual trajectories throughout the simulation that helped to gain further insight into the DOR complexes (Table S5, Figure 5, Figures S10-S12). Expecting the agonist bound DOR to undergo conformation change, it can be understood that the majority of the individual agonist trajectories contained more than one populated cluster conformation. In contrast, the individual antagonist trajectories each generated only one populated cluster conformation, as they were not expected to have changed from their initial inactive conformation.
After calculating the RMSD of different portions of the receptor when superimposed with the crystal structures, we were able to identify three independent activating events of the ADL5859 system that supports the DOR folding toward the active conformation, with the most notable hallmark as the RMSD of TM6. The first event occurred in the first trajectory where cluster 1 (52%) and cluster 3 (17%) were structurally similar to the active crystal complex with favorable RMSD values (Cluster 1 RMSD TM6 Inactive: 4.1 Å Active: 3.4 Å Cluster 3 RMSD TM6 Inactive: 3.3 Å Active: 2.9 Å). The second event happened in the fourth trajectory where cluster 5 (1%) showed agreement with the active crystal complex with favorable RMSD values (Cluster 5 RMSD TM6 Inactive: 5.8 Å Active: 2.8 Å). The third event was in the sixth trajectory's most abundant cluster (76%) that was structurally similar to the active crystal structure with favorable RMSD values (Cluster 1 RMSD TM6 Inactive: 3.9 Å Active: 3.2 Å). The percentage of abundance of each event indicates >70% folding into the active state with the 1% of abundance from trajectory 4 cluster 5 as a partial folding event. The agonist bound clusters shared similar interacting protein-ligand residues as the active crystal structure as well. In the three individual naltrindole trajectories, each abundant cluster showed structural agreement with the inactive crystal structure (4N6H) with lower RMSD values against 4N6H and higher values against 6PT3, further inferring its similarity to the inactive conformation. This is further validating that the activation related changes are unique to our agonist bound complex.
Superimposing the abundant structures to the crystal structures and calculating the RMSD of different protein regions also opened speculation to other populated clusters from the individual agonist trajectories as intermediate conformations in the protein folding toward the active state. Trajectory 3 clusters 1 (90%) and 2 (11%); trajectory 4 clusters 1 (43%), 3 (20%), and 4 (5%); and trajectory 5 cluster 1 (56%) show some to great agreement to the active crystal structure, especially in regions TM5-6, while showing disagreement to the inactive crystal structure in these same regions. These speculated intermediate active state folding structures add further support that our agonist-ADL5859 system has undergone protein folding toward the active state as the trajectories from which the structures came (trj3, trj4, trj5), in addition to the respective trajectories of the identified active folding structures from the clustering RMSD calculation (trj1, trj4, trj6), are consistent with the individual agonist trajectories determined as adopting toward the active conformation from the RMSD time series of the whole simulation (trj1, trj3, trj4, trj5, trj6) ( Figures S8 and S9). This analysis has not only offered support to the previous RMSD analysis, but has also added more insight into the speculated activation process towards the activated DOR.
Comparing the two representative structures from the MD complexes revealed a clear conformational difference in the receptors in the transmembranes as well as the intra-and extracellular loops and the ligand binding pose. These differences in the protein structures are consistent with the differences between the active and inactive crystal structures. The binding pose of ADL5859 is lower than that of naltrindole, which may be a result of the great conformational change occurring that could be pulling the ligand deeper into the binding pocket. The MD simulations were successful in producing different conformations in the DOR when bound to an agonist and antagonist. Our simulation data supports the conformation selection mechanism for DOR activation.
The RMSF of the receptor (Figures 7 and S20-S21) when bound to agonist-ADL5859 indicated the complex had greater flexibility in the loops and terminals when compared to the antagonist-naltrindole system which is consistent with not only the pharmacological actions of the two ligands (agonist/antagonist), but also our hypothesis that our agonist system adopted toward the active conformation. The ligand RMSF also supports our hypothesis due to ADL5859 itself showing greater flexibility than naltrindole. This is consistent with what is known, that an antagonist (naltrindole) typically binds tighter to the receptor than an agonist (ADL5859). Naltrindole is also a more rigid ligand compared to ADL5859, with more aromatic rings and less rotatable bonds, which tells enough in itself about the predicted flexibility of the ligands that was produced in the ligand RMSF. The RMSF analysis is another indication supporting that our agonist system underwent a conformational change.
Residue D95 2.50 is believed to be an important residue in the agonist-induced cyclic adenosine monophosphate (cAMP) response that occurs in activated GPCRs (White et al., 2018). In timeline plots, we see that in traj1 and traj6 decrease in the distance between ADL5859 and D95 correlates with DOR activation (Figure 8). In contrast, the distance between D95 and naltrindole has no change throughout the simulation, inferring that the helices have not shifted, remaining in the inactive state.
Molecular switches are hallmarks in the classification of the activation state of class A GPCRs (Trzaskowski et al., 2012). Previous studies have identified the lack of a salt bridge in inactive opioid receptors which have a hydrophobic Leu 6.30 in place of the usual acidic Glu 6.30 (Fenalti et al., 2014). It was later realized hydrogen bonding between R146 3.50 and T260 6.34 residues may stabilize the receptor in the inactive state similar to the salt bridge in the Ionic Lock Switch (Katritch et al., 2013). This interaction was analyzed as the Ionic Lock Switch along with the usual Transmission and Tyrosine Toggle Switches. This hydrogen bonding was broken in multiple trajectories with a distances of 4.2-9.3 Å. A very recent study by Claff and coworkers reported on the distance of the transmission switch on a mutant active state DOR with a large outward movement on helix 6 with a distance of approximately 9.4 to 11.2 Å (Claff et al., 2019), which is consistent with our data that showed distances of 9.2-13.9 Å. The sidechain of W274 of the transmission switch also showed slight rotameric change, indicating the receptor is in the active state. Any slight change in the orientation of side chain of W274 can have direct effect on the residues of TM6 and other regions of DOR. This same study as previously mentioned also reported an inward shift of helix 7 at a distance of 3.9 Å, however, they speculate this greater shift may be attributed to the three mutations they made in the sodium-binding pocket (N90 2.45 S, D95 2.50 G, N131 3.35 S) (Claff et al., 2019). In contrast, this shift was not observed in the trajectories of the agonist-ADL5859 bound DOR. Tyrosinetoggle switch did not show significant change in ADL5859-DOR complex. As for the Naltrindole complex, the three molecular switches remain intact, which obstruct the outward movement of TM6 thus, blocking G protein signaling.
In addition, we have used multiple atomic and center of mass distances to characterize the conformation differences between the active and inactive states of the DOR. These differences could contribute the activation. Indeed, the seven residue pairs were identified that are associated with activation.
The normal mode analysis ( Figures S33 and S34) shows the top 10 low vibrational modes identified from the principal component analysis. Previous knowledge tells us that the motion and conformation of transmembrane six is relative to the state of the receptor (i.e. active, inactive). While the overall motion is significant, this area can be focused in on for clarity. We conducted this analysis on trajectory 1 and 6 from our agonist system and all three antagonist trajectories. Our analysis identified the ADL5859 system having a higher degree of fluctuation and a more distinct motion when compared to naltrindole. Altogether, the two agonist trajectories had seven conserved vibrational modes in which they had similar motion (mode 1, 2, 3, 4, 7, 8, 9). The conserved vibrational modes where motion was favored between the two trajectories were modes 1, 2, 4, and 8 where the ends of TM5 and TM6 swing outward, opening the intracellular portion of the receptor for G protein signaling that is crucial in GPCR activation. The motion of the naltrindole system in the majority of the modes for all three trajectories of the same region can be observed as moving back towards TM3 and TM4 or downward, still in a closed inactive position. This further provided another level of insight into the essential dynamics of these complexes in the lowest vibrational states. We can conclude from the normal mode analysis that the overall dynamics of the two complexes are distinct from one another, which further supports that the motion differences of these complexes are due to their dynamic interactions as well as their pharmacological actions. Our analyses including RMSD, RMSF, molecular switches, measuring the center of mass, and normal mode analysis help to explain the activity differences on the receptor that cannot be explained from the crystal antagonist alone.
The proposed activation pathway of the receptor started from the allosteric signal from the agonist ligand to the intracellular domain that interact with the effector G-protein has been summarized in Figure 10. When the agonist is bound to the DOR, it forms strong interactions with W274 and D95 in layer 1-2, further leads to break important interaction pairs from layer2 by layer4 such as V267-L139, V263-M142 and others, and finally causing an outward shift of the intracellular domain of TM6. This is consistent with what is known about class A GPCRs that the intracellular portion of this helices should move outward when activating to allow for G-protein signaling (Trzaskowski et al., 2012). In this same respect, TM3 shifts away from TM6, breaking any hydrogen bonding that stabilizes the inactive conformation. While this occurs, TM7 also moves away from both TM3 and TM6. During all of this movement, the extracellular portion of TM1 folds away from the ligand and the binding pocket while TM2 shifts closer to the ligand. The experimental validation of the proposed activation mechanism is required, mutagenesis studies could be a good tool for such purpose. The key residues in the activation pathway identified from this study provides a short list for future experimental mutagenesis study.

Conclusion
The crystal structures of the DOR in complex with an antagonist and agonist have been recently solved, however, the active state contained mutations in the receptor. Using a different agonist that has no crystal structure and the crystal structure of the antagonist naltrindole has provided dynamic insight into the activation mechanism of the DOR. Comparing the trajectories and protein-ligand interactions allowed us to gain a better understanding of the activity and conformational changes induced from each ligand. While the protein-ligand interactions and secondary structures showed subtle variations, the protein RMSD in reference to both crystal structures (4N6H & 6PT3) for both systems were where the major differences were observed; with the ADL5859 complex showing to adopt toward the active conformation and the naltrindole complex showing to remain in the inactive conformation. In addition to this, the protein and ligand RMSFs for both ligand systems showed great differences as well. ADL5859 displays a trend of higher RMSF values, indicating it to have higher fluctuation than naltrindole. Observing the structures of the two ligands can explain this with naltrindole consisting of many rings, making it a more rigid structure with less flexibility than ADL5859. Interaction of ADL5859 with transmission switch, particularly W274 appeared important for the DOR activation. By examining the conventional molecular switches and additional pairs of interactions, the key residues such as W274, D95, V267, L139, V263, M142, I88, T260, R146, R258 and others in the activation pathway were identified, providing a short list for future experimental mutagenesis study. The differences in dynamics between the two complex systems was explored and made more clear by using normal mode analysis. All of our results and key findings suggest that the conformation selection mechanism is the activation mechanism of the DOR by agonist ADL5859. Further experiments such as mutagenesis study are required to validate our proposed activation mechanism.

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