Probing the Activation Mechanisms of Agonist DPI-287 to Delta-Opioid Receptor and Novel Agonists Using Ensemble-Based Virtual Screening with Molecular Dynamics Simulations

Pain drugs targeting mu-opioid receptors face major addiction problems that have caused an epidemic. The delta-opioid receptor (DOR) has shown to not cause addictive effects when bound to an agonist. While the active conformation of the DOR in complex with agonist DPI-287 has been recently solved, there are still no FDA-approved agonists targeting it, providing the opportunity for structure-based virtual screening. In this study, the conformational plasticity of the DOR was probed using molecular dynamics (MD) simulations, identifying two representative conformations from clustering analysis. The two MD conformations as well as the crystal conformation of DOR were used to screen novel compounds from the ZINC database (17 million compounds), in which 69 drugs were picked as potential compounds based on their docking scores. Notably, 37 out of the 69 compounds were obtained from the simulated conformations. The binding stability of the 69 compounds was further investigated using MD simulations. Based on the MM-GBSA binding energy and the predicted drug properties, eight compounds were chosen as the most favorable, six of which were from the simulated conformations. Using a dynamic network model, the communication between the crystal agonist and the top eight molecules with the receptor was analyzed to confirm if these novel compounds share a similar activation mechanism to the crystal ligand. Encouragingly, docking of these eight compounds to the other two opioid receptors (kappa and mu) suggests their good selectivity toward DOR.


■ INTRODUCTION
Pain signal processing, or pain perception, occurs when the brain is alerted by electrical signals transmitted by nociceptors due to tissue damage.These electrical signals are transmitted to the central nervous system to alert the brain about potential harm. 1 The key membrane receptors involving pain signal processing are nociceptors, sensory neurons, glial cells, and postsynaptic neurons.As previously mentioned, nociceptors are responsible for detecting potential damaging stimuli at the skin caused by environmental factors, such as chemicals or heat. 1 Sensory neurons, such as sodium and potassium channels, are responsible for transmitting pain signals from the wounded area to the central nervous system. 1 Glial cells, such as toll-like receptors (TLRs), are responsible for releasing the inflammatory mediators.Postsynaptic neurons, which are the opioid receptors, provide analgesic effects when activated.
In pain perception for opioids, the neurotransmitter that is linked to the pain sensation in the central nervous system (CNS) and peripheral nervous system (PNS) is called substance P (SP).SP is a neuropeptide whose role is to act on its receptor, known as the neurokinin-1 (NK1) receptor, to transmit pain signals in the CNS.Nociceptors are responsible for releasing SP following tissue damage. 2 When SP activates the NK1 receptor, it releases inflammatory mediators, such as histamine and cytokines, from the peripheral nerves. 3This process is known as neurogenic inflammation, which enhances pain signals and promotes tissue inflammation. 4During pain sensation, the brain signals the release of endorphins to relieve pain and other opioid-associated effects, such as respiratory depression.Endorphins are endogenous opioid peptides that the body produces to bind to opioid receptors for regulating pain perception. 5This also leads to the activation of internal opioid receptors caused by external opioids.As a result, it leads to pain relief, euphoria, and respiratory depression.There is no clinical evidence to support that the NK1 receptor antagonist is a good pain-killing drug target in pain management.
Some major pain medicines used to treat pain are opioids, anticonvulsants, and topical analgesics, such as gels.Anticonvulsants are used in neuropathic pain management; they bind to specific subunits of voltage-gated calcium channels to regulate neurotransmitter release and reduce neuronal activity. 6n topical analgesics, sodium channels are blocked to reduce pain signals. 7Nonetheless, opioid drugs are essential in pain management.
While pharmaceuticals have made great strides over the years in many areas, the field of pain management is still lacking.There are great medications that target pain and are effective; however, there has been a rise in the use of opioids over the years that has led to the opioid epidemic in Northern America and other parts of the world.The use of opioids (legally and illegally) has risen between 10 and 14 times in the last 20 years, with more than 42,000 deaths in 2016 in the USA occurring from opioid overdoses alone. 14,15Most prescribed opioids target the muopioid receptor (MOR) in Table 1, which are located in the reward areas of the brain. 14,16,17When opioid agonists bind to and activate these receptors, it causes euphoria, which can lead to addiction and severe withdrawal symptoms after repeated use.The need for better pain management without addictive properties is thus pressing.
−25 Due to the physiological symptoms that occur when opioid agonists bind to the MOR, this makes the DOR an attractive target to further study to potentially help alleviate the opioid epidemic in the world.Like the DOR, the kappa-opioid receptor (KOR) has also shown to not have addictive properties and does not induce respiratory depression. 11When this receptor is activated by its endogenous ligand dynorphin, it produces analgesic, physiological, and behavioral effects (Table 1). 12A study has shown that the selective KOR agonists, such as U50, can be considered as a nonaddictive alternative without the analgesic effects. 13here are multiple studies that have reported convulsions in various animal models with the use of DOR agonists.After systemic administration of a DOR agonist, mice displayed convulsive effects. 8Using rats, tolerance rapidly developed to convulsive and locomotor-stimulating effects of a selective DOR agonist but did not display tolerance to the antidepressant-like effects. 9When using rhesus monkeys, only one out of the four monkeys had convulsions; however, this same monkey did not display convulsive activity when given a smaller dose weeks later or even the same dose 1 year later. 10The difference in convulsions in species could indicate that these are speciesdependent effects, and further studies are required to resolve this issue.As mentioned previously, the DOR is distributed in different areas of the spinal cord in rodents versus primates.In rodents, the delta receptor is found dispersed in the spinal cord whereas it is limited to the superficial laminae of the spinal dorsal horn in humans and non-human primates.
Although DOR agonists have displayed convulsive effects in mice, and rhesus monkeys after administration, the DOR is still an attractive target for treating chronic pain, anti-depression, and ischemic preconditioning (Table 1).DPI-287 is an experimental opioid agonist found to induce less convulsions than most drugs in this family and is more selective toward the DOR than MOR or KOR in the rat forced swim test. 58,26The DOR has a possible site where morphine tolerance and dependence can be regulated. 61Further analysis showed that DPI-287 had weaker docking scores at MOR and KOR and stronger docking scores at DOR, demonstrating a 10-fold selectivity increase for DOR.With DPI-287 being more selective toward the DOR and having reduced risks for convulsions, it makes the DOR an attractive target.Despite this, no FDAapproved drugs currently exist for targeting the DOR; thereby, discovering potential lead compounds for DOR treatment is of utmost importance.
The integration of high-throughput virtual screening (HTVS) in drug discovery is very useful in screening for thousands of molecules that bind favorably to a molecular target and for identifying toxic or unfavorable pharmacodynamic and pharmacokinetic properties of these compounds. 27Molecular dynamics (MD) simulations are important in drug discovery due to the significant kinetic and mechanistic information it provides. 64Structure-based virtual screening (SBVS) is an HTVS approach that predicts the interactions between ligands and proteins as a complex, ranking them by their affinity to the receptor.The top hit compounds are then selected based on the desired parameters and are then optimized to undergo preclinical and clinical trials.Other computational methods such as molecular modeling are used in the HTVS approach.The process speeds up the drug discovery development by analyzing the interactions of multiple molecules in a shorter period of time, which can look into interactions before the drug is even synthesized.Hence, SBVS is a good technique due to its low cost, faster result time, and good results achieved.
The flexibility of receptors is a challenge researchers face as binding sites usually consist of 10−20 amino acid side chains that have multiple rotatable conformations, which are larger than the rotatable torsions of a ligand. 28The movements of the protein backbone can make this even more complex by affecting the conformations of multiple side chains.Using an ensemblebased receptor technique combats this issue in HTVS and MD simulations by sampling the degrees of freedom instead of traditional techniques using one receptor conformation with a flexible ligand.It has been found in cases to also improve docking scores.With previous studies, it has been suggested that ensembles generated from MD simulations have been closely similar in replicating the dynamics of proteins in NMR experiments. 29,30It is better to use a few specifically selected conformations as using too many could give false results.In the case of virtual screening, using the top 10% of a compound library subset is more efficient with this approach. 28Previous studies using ensemble-based virtual screening have been successful in screening ligands against various drug targets and provide insight for future drug designs. 31,32,65The integration of this ensemble-based technique helps to have a better understanding of the structural dynamics of a receptor and have a better understanding of ligand−receptor interactions, which aids in discovering novel ligand binding modes, and helps to develop better therapeutic molecules. 33ynamic network analysis was used to interpret and understand atomic information and structural changes among different regions of the DOR. 63,66Unweighted and weighted network models were calculated to decipher the allosteric signal transmission pathway, and this comparison showed their connection to be in good agreement.With MD simulations and network models, this allowed for identification and comparison with interactions of signaling pathways that occur in a system. 62Ensemble-based screening was also used in this study, and it is a cost-effective method that has the potential to provide breakthrough discoveries, by providing accurate estimates of free binding energies.To the best of our knowledge, this is the first study to use ensemble-based screening to discover potential agonists to target the DOR.
In the present study, MD simulations were used to probe the active conformation of the DOR starting with the active crystal conformation (PDB ID: 6PT3) and the crystal agonist DPI-287 in Figure 1.Using the ensemble-based method, two representative conformations were identified from the clustering and principal component analysis based on the MD simulations.These two conformations and the crystal conformation were used to screen 17 million compounds from the ZINC database.As a result, 69 compounds were identified based on their docking scores.These 69 protein−compound complexes underwent MD simulations to assess their stability.From this, eight compounds were identified that showed significantly improved MM-GBSA binding free energy scores with high blood−brain barrier (BBB) permeability and high gastrointestinal (GI) absorption.This study helps to identify potential compounds to be further tested that will aid in antinociception without addictive or convulsive properties for the DOR.

■ METHODS
Using the ZINC15 "drug-like" library, which contains 17 million entries including FDA-approved drugs (∼1615 entries), a virtual screening workflow (VSW) was developed to identify lead agonists to the DOR (Figure 2) The ZINC library defines "druglike" using a widely used Lipinski's Rule of 5: molecular weight (≤500), number of hydrogen bond donors (≤5 H-bonds from the sum of NH and OH), number of hydrogen bond acceptors (≤10 H-bonds from the sum of N and O), and the octanol− water partition coefficient LogP (≤5). 34,35The VSW comprises 10 steps including molecular docking, drug property prediction, MD simulations, and post MD-simulation analysis.Inputting the prepared protein structure and ligand library is the first step of the VSW.In steps 2−5, the compounds were then filtered by docking scores and multiple Glide docking score functions that have increasing accuracy (Glide HTVS, SP, and XP).In the next step, a ligand similarity analysis was performed to identify different molecular scaffolds.In step 7, the ligands removed were based on if they had a worse Glide XP score than the reference compound DPI-287 (PDB ID: 6PT3) and/or if they had more than one red flag in drug property (number of stars, from QikProp).QikProp is a Maestro software module that evaluates the drug-like characteristics of a compound by comparing its molecular properties to those of established drugs.Unlike Lipinski's rule, which uses only four properties, QikProp examines 24 molecular properties.The drug likeness of a specific molecule is forecasted by applying a Gaussian probability distribution function to each of these 24 properties, drawing from data on FDA-approved drugs.For each property, the cutoff value is defined as the value that 95% of FDAapproved drugs adhere to.A compound is awarded a violation star (*) across all of its properties, and the total number of stars is then calculated.As such, fewer stars (0−1 star as per this study) suggest fewer violations and toxicity, indicating a compound is more drug-like.A list of the 24 properties analyzed by QikProp can be found in the supporting document (Table S1). 36he top compounds were manually selected from the remaining compounds by maximizing the number of molecular scaffolds (i.e., different ligand cluster IDs).The 200 ns MD simulations were carried out in step 8, followed by the post-MD simulation analyses in step 9, including MM-GBSA binding free energy calculation, simulation interaction diagram analysis, and protein conformation clustering analysis.In the last step, the ADMET (absorption, distribution, metabolism, excretion, and toxicity) prediction was used to examine the human oral bioavailability of potential drug candidates.From this, the compounds with better MM-GBSA binding free energy compared to the reference compound were selected and presented in the main text.The 10 steps are presented in detail in the following six modules.
1. Preparation of the Protein and Ligand Library.The crystal structure (PDB ID: 6PT3) of the active conformation of the DOR was prepared and preprocessed using Maestro's Protein Preparation Wizard. 37The preprocessed protein's charge state was optimized at pH 7.Then, a restrained minimization was performed to relax the protein structure using an OPLS3 force field. 38Maestro Elements was used to prepare the 3D structures of DPI-287 and the ZINC compounds and FDA-approved drugs (2466 entries).The 3D structure of DPI-287 was extracted from the crystal structure (PDB ID: 6PT3), and the ZINC compounds and FDA-approved drugs were downloaded from the ZINC15 database (https://zinc15.docking.org/)and the Drug Bank, respectively.In order to generate each ligand's ionization/tautomeric states at pH 7, Maestro's Epik tool was used, based on the Hammet and Taft methodologies for increased accuracy. 37The lowest ionization/ tautomeric state was then chosen.Afterward, the ligand geometry was optimized using quantum mechanics in Maestro's Jaguar tool.
2. Filtering and Docking.The prepared protein and ligand structures were merged into a complex to then be ran through the Schrodinger Virtual Screening Interface using prefilters through Lipinski's Rule and filtered with ADMET risk parameter assessments through QikProp.The receptor grid files were generated from the prepared receptors, in which the centroid of the crystal ligand, DPI-287, was used to specify the active site.The prepared ligands were docked into their corresponding generated grids using Glide XP scoring with default procedures and parameters. 38In detail, the receptor grid required for the docking process was generated using a van der Waals scaling factor of 1 and a partial charge cutoff of 0.25.Docking was performed using a ligand-centered grid and an OPLS3 force field.Glide XP Dock was used perform a comprehensive systematic search for the best receptor conformations and orientations to fit the ligand.The docked poses were compared to the active crystal complex (PDB ID: 6PT3) with the agonist to verify if the docked ligand poses were reasonable.All ligands were bound within the binding pocket with DPI-287 binding similar to the crystal ligand, 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.The docking results comprised 69 top ZINC compounds and top 10 FDA-approved compounds.Notably, the 69 top ZINC compounds exhibited much higher docking scores than the reference ligand (PDB ID: 6PT3) and the top 10 FDA-approved compounds, indicating that these ZINC compounds all had high affinity for the receptor.Additional analyses were thus only done for the top 69 ZINC compounds.
3. Ligand Similarity Clustering.The Canvas program was used for determining ligand similarity clustering.Canvas uses pharmacophore fingerprinting and hierarchical clustering to further filter out the top 69 ZINC compounds.Pharmacophore fingerprinting identifies similar groups of compounds to match the crystal structure. 39Hierarchical clustering forms cluster groups of similar compounds based on their docking score, binding affinity, drug properties, and ligand similarities. 40,41ncouragingly, the top 69 ZINC compounds were then used for the MD simulations.
4. MD Simulation.4.1.MD Simulation System Setup.The 69 prepared receptor−ligand complexes from Glide XP docking were used as input files to construct MD simulation systems using Desmond System Builder. 38Each complex was inserted in a phosphatidylcholine (POPC) lipid membrane, 42 solvated in an SPC water box of orthorhombic shape with a 6 Å water buffer, 43 and Na + and Cl − ions neutralized the intrinsic system charge and established a 0.15 M NaCl salt concentration using the OPLS3 force field.

Relaxation and Production Runs.
A default relaxation protocol for membrane proteins in the Desmond module was used to relax the system.The relaxation protocol consisted of eight steps.The first step minimizes the solute heavy atoms with restraints.The second step minimizes the solute heavy atoms without restraints.The third step performs a heat transfer simulation from 0 to 300 K, which led to a H 2 O barrier and gradual restraining.The fourth step performs a simulation under the NPT ensemble (constant pressure (1 bar) and temperature (300 K)) with a H 2 O barrier and heavy atoms becoming restrained.The fifth step performs a simulation under the NPT ensemble to equilibrate the solvent and lipid components.The sixth step performs protein heavy atoms annealing from 10.0 to 2.0 kcal/mol by performing a simulation under the NPT ensemble.The seventh step was to restrain protein Cα atoms at 2 kcal/mol by simulation under the NPT ensemble.The eighth step was to execute simulation for 1.5 ns under the NPT ensemble with no restraints.
After the relaxation protocol, each equilibrated MD system was put through a 200 ns production run under the NPT ensemble using the default protocol.During this process, the Nose−Hoover chain coupling scheme 44 was used to control the temperature at a coupling constant of 1.0 ps.The Martyna− Tuckerman−Klein chain coupling scheme 44 was used to control the pressure at a coupling constant of 2.0 ps.To restrict all bonds connecting hydrogen atoms, M-SHAKE 45 was applied, enabling a 2.0 fs time step in the simulation.The k-space Gaussian split Ewald method 46 was used to handle long-range electrostatic interactions under periodic boundary conditions.During this step, the charge grid spacing was set to about ∼1.0 Å while the direct sum tolerance was 10 −9 .The short-range non-bonded interactions had a cutoff distance of 9 Å.Meanwhile, the longrange van der Waals interactions were established on a uniform density approximation.An r-RESPA integrator 47 was used to reduce the computation by calculating the non-bonded forces.The short-range forces were updated at each step, while the long-range forces were updated every three steps.The trajectories at 50.0 ps intervals were saved for further analysis.
5. Post-Simulation Analysis.5.1.Simulation Interaction Diagram (SID) Analysis.The Desmond SID tool was used to generate protein/ligand root mean square deviation (RMSD), protein−ligand contacts, secondary structure changes, and protein/ligand root mean square fluctuation (RMSF) measures.The protein Cα and ligand heavy atom RMSD plots were analyzed to check the convergence of MD simulations (i.e., steady-state equilibrium).
5.2.Trajectory Clustering Analysis.For each complex system, the Desmond trajectory clustering tool 48 was used to combine complex structures from the last 100 ns of each simulation.The backbone RMSD matrix was used as a structural similarity metric, while the hierarchical clustering with average linkage was used as the clustering method.The merging distance cutoff was 2 Å.The centroid structure was represented as the structural family, where the centroid structure represents the largest number of neighbors in the structural family.

Binding Energy Calculations and Decompositions.
The surface-area-based generalized Born model 49,50 was used to calculate the ligand-binding affinities on the frames, in the last 50 ns of each MD simulation, with an implicit membrane solvation model (VSGB 2.0). 51Slab-shaped regions with a low dielectric constant between 1 and 4 were excluded from the implicit membrane and were assigned with the solvent (water) dielectric constant of 80.The MM-GBSA calculation used an OPLS3 force field and the default Prime procedure. 38The OPLS3 force field employs a CM1A-BCC-based charge model based on a combination of Cramer−Truhlar CM1A charges 52 with an extensive parameterization of bond charge correction (BCC) terms.This process begins with minimizing the receptor only, then the ligand only, and then the receptor−ligand complex.Using eq 1, the MM-GBSA binding free energy for each system was calculated from three separate simulations: ligand-only, receptor-only, and the receptor−ligand complex.Equation 2 contains four components: van der Waals interaction energy (VDW), hydrophobic interaction energy (SUR), electrostatic interaction (GBELE), and the change of the conformation energy for the receptor and ligand that were calculated based on eqs 3 and 4. (1) The MM-GBSA scoring function lacks the solute conformational entropy contribution, which causes higher negative values when compared to actual values.It is essential to rank a drug's ability to target a receptor when it is used to rank different drugs targeting receptors with comparable entropy values. 53−59 6. ADMET Prediction.ADMET properties were predicted for the best ZINC compounds and were performed on the SwissADME web server (http://www.swissadme.ch/).The SwissADME server was developed by the Swiss Institute of Bioinformatics and is used to provide physiochemical descriptors, ADMET parameters, pharmacokinetic properties, and drug-like small molecules to support drug discovery. 60In order to receive each compounds ADMET properties, their respective SMILE codes were inserted into the web server as inputs.

Normal Mode Analysis.
The Normal Mode Wizard plugin in VMD 61 with default settings was used to generate a main component analysis of the top 10 normal modes for the crystal complex (PDB ID: 6PT3).
8. Dynamical Network Model.A dynamic network model, defined as a set of nodes connected by edges, 58,62−65 was generated using the individual trajectories of each system using the NetworkView plugin in VMD. 66We first generated a contact map for each system of the top compounds and the crystal complex, which added an edge between nodes, whose heavy atoms interacted within a cutoff of 4.5 Å for at least 75% of the MD simulation time.The edge distance was derived from pairwise correlations 64 in the contact map using the program Carma, 67 which defines the probability of information transfer across a given edge using the following equation: (5)   The edges in the dynamic network model are weighted (w ij ) between two nodes i and j, which uses the following calculation: The weight of the edge is correlated with the probability for information to transfer across the edge between two nodes.Because of this, a thicker edge is characterized as a higher probability of information transfer.The network for each system was further grouped into communities, or subnetworks based on groups of nodes with more frequent and stronger connection to each other, by applying the Girvan−Newman algorithm to the original network. 68The critical nodes that connect communities to each other were identified as well.Optimal communication paths were generated between the ligand node and the molecular switch residue number using the data from the molecular switches.9. Selectivity Analysis on Different Opioid Receptors.To determine the relative selectivity of the top eight ZINC compounds to DOR, Glide XP docking of the top eight ZINC compounds was also performed on the KOR and MOR structures.The MOR is taken from our previous studies, 69,70 which is based on the human MOR homology model built from the crystal structure of mouse MOR in complex with the agonist BU72 (PDB ID: 5C1M).The KOR is taken from the PDB bank (PDB ID: 6B73) and is based on the crystal structure of a nanobody-stabilized active state of KOR.The Glide XP docking scores of the eight ZINC compounds were compared between those of the DOR crystal structure and of the KOR crystal structure and MOR homology model structure.Furthermore, the crystal ligands DPI-287/DOR, MP1104/KOR, and BU27/ MOR were also docked to each of the three receptors to serve as reference values.

Stability of Crystal Conformation of DOR Maintained
during MD Simulation.The crystal complex of the active conformation of the human DOR (PDB ID: 6PT3), with crystal agonist DPI-287, was used in the experiment to serve as the control.DPI-287 was first docked back into the crystal conformation and resulted in a similar binding pose with the crystal ligand pose (RMSD = 0.47 Å) and a docking score of −8.6 kcal/mol.The low RMSD value validates the docking protocol (Figure S6, row 1).A 1000 ns MD simulation production run and post-MD simulation analyses were then performed to check the stability of the prepared crystal structure complex.Indeed, the complex was stable throughout the entire trajectory (Figure 3).
To check if the simulation system reached convergence, the protein and ligand RMSD was calculated over the whole 1000 ns trajectory (Figure 3A).The last 200 ns shows that the protein and ligand system components converged with relatively flat plots.The average RMSD values of the protein (2.9 ± 0.1 Å) and ligand RMSD (2.5 ± 0.1 Å) are relatively small, suggesting the crystal complex was stable throughout the trajectory (Table 2).
To determine any changes in the protein secondary structure, secondary structure element (SSE) plots were generated for each protein residue.The SSE plot summarizes the structure element distribution by residue position throughout the protein.
The three categories are alpha-helices, beta-strands, and random coils.Alpha-helices are mainly made up of hydrophobic residues that are located in the core of the protein and are depicted by the orange sections.Beta-strands, however, contain both hydrophobic and polar amino acids, which are depicted by blue sections.The random coil is not one specific shape of a polymer conformation but a distribution of statistics of all chains depicted by the white spaces in the plots.SSE indicates that the helices were maintained during the simulation (Figure 3B).
RMSF was then calculated to determine fluctuations in each protein residue, intra-and extracellular loops, and N-and Ctermini, which are usually the most flexible regions of the receptor.While some protein regions exhibited some increased RMSF spikes, the average RMSF was relatively small (<2 Å) and indicates low global fluctuations (Figure 3C).Finally, the 2D ligand interaction diagram was generated to determine the specific interactions between the crystal ligand and the protein active site, which revealed mostly hydrophobic contacts and some hydrogen bonding (Figure 3D,E).Overall, the system was shown to be stable and mimic the crystal structure with a stable binding pose.
Crystal Complex Produces Other Conformations to Use for HTVS.Clustering analysis was done after MD simulation to identify the populated conformations for the trajectory.Each cluster conformation contains a percentage of abundance based on the clustering algorithm in which a cutoff of 2% was used.From this, there were two abundant clusters (75 and 24%, respectively) produced from the crystal conformation simulation and were compared to the crystal complex (Figure 4).The two clusters slightly differ in conformation and the ligand binding pose from the crystal conformation.In a more precise view, the binding pocket of each structure was compared to the crystal pose (Figure 5).In this view, differences can be seen in the receptor itself and the side chains that have adopted different rotamer states.Specific residues where the side chains differ the most from the crystal in both clusters are N90, D95, D128, N131, N310, and Y318.
Normal mode analysis on the crystal MD simulation trajectory revealed the top two low-frequency vibrational modes (Figure 6).This further validates the result of two cluster conformations that differ from the crystal conformation    and therefore can be used as additional conformations to use for HTVS.Unweighted and weighted dynamic network models of the DPI-287/DOR system were calculated as described in the Methods section to decipher the allosteric signal transmission pathway.The unweighted network model shows that their connections are in good agreement with each other.Quantification of the correlation between the nodes in the weighted network model reveals the areas of the receptor that are in higher correlation to each other.The system appears to have higher correlations between edges TM5 and TM6.
A community network model was generated using the weighted network model, which grouped residues together that interact more frequently and stronger than to residues in other communities (Figure 7).There were 38 critical nodes identified for the crystal system from the community analysis.These critical nodes were involved in signal transduction between different parts of the receptor throughout the simulation; therefore, the critical residue information was then cross referenced with experimentally reported mutagenesis data available on the G protein-coupled receptor databank (GPCRdb) to see if the residues were involved in the physical signal transduction.The DPI-287 system had nine critical residues that were also naturally occurring mutations (D95, Y129, M132, M142, R160, M186, A269, I282, N310).It also had three critical residues that were mutations in vitro (D95, W274, S312).Optimal paths generated for the DPI-287/DOR system give insight into the molecular signal transduction pathways involving the ligand.From the weighted network models, the shortest pathways able to pass a signal from the ligand to the site of the molecular switch (Tyrosine Toggle Switch: Y318) and the intracellular end of TM6 (Transmission Switch: S255) were calculated as the optimal paths.DPI-287 has a direct optimal path for the Transmission Switch (CWXP) through TM6 and another direct path for the Tyrosine Toggle Switch (NPXXY) through TM7.
Top Eight ZINC Compounds Identified from MD Simulations, MM-GBSA, and SwissADME Property Prediction.HTVS was ran on all three conformations of the active DOR (crystal conformation, cluster 1, cluster 2) for FDAapproved 2466 drugs 32 and compounds from the ZINC15 database, as described in the Methods section.Using the docking score of the crystal ligand DPI-287 as a cutoff for favorable binding (−8.6 kcal/mol), the number of stars (0−1 stars) for drug-likeness, and the cluster ID in picking diverse chemical scaffolds, a total of 69 ZINC compounds were chosen.The top 10 FDA-approved drugs were also identified from their docking scores.However, comparing the docking scores revealed that the top 69 ZINC compounds had significantly better docking scores than the FDA-approved drugs, therefore ruling out the further investigation of the FDA-approved drugs (Figure S2, Table 2, and Table S2).From the 69 ZINC compounds, 32 were docked to the crystal receptor (Figure S3), 11 were targeted to the first abundant cluster receptor (Figure S4), and 26 were targeted to the second abundant cluster receptor (Figure S5).Each docked ZINC compound showed good binding pose agreement with the crystal ligand (Figures S6−S8 MM-GBSA was calculated for each of the 69 ZINC compounds to their respective receptor conformation, including the crystal ligand DPI-287, to determine the binding affinity for the receptor, with more negative values indicating better binding (Table S2).DPI-287 had an MM-GBSA of −90.2 kcal/mol, which was used as the cutoff and reference.From this, the top 69 ZINC compounds were filtered down to 22 compounds, based on their favorable MM-GBSA binding energies (−91.5 ± 7.8 to −134.9 ± 11.3 kcal/mol).
To further filter down the top ZINC compounds, SwissADME properties were then calculated for the top 69 ZINC compounds (Table 3 and Table S3).SwissADME properties include but are not limited to GI absorption and BBB permeability, which were the major determinates.Also, similar to the crystal ligand, having no alerts (PAINS, BRENKS) was attractive as well due to having a low chance of false positives from occurring.Considering these desired properties, eight ZINC compounds were chosen as the top compounds (Tables 2  and 3), which also exhibit more favorable MM-GBSA binding free energies (Table 2).Two of the ZINC compounds were targeting the crystal conformation, four were targeting the first abundant cluster, and two were targeting the second abundant cluster.From here on, the top eight ZINC compounds are the focus of additional analysis.
Top Eight ZINC Compounds Assume Steady-State Equilibrium.Protein and ligand RMSDs were calculated over each trajectory for the top eight ZINC compounds to check for convergence (Figure 8).RMSD plots for the remaining 61 ZINC compounds are in the supporting document (Figure S9).All eight MD simulations achieved steady-state equilibrium during the last 50 ns of simulation time, with the lower protein RMSD values indicating they remained stable throughout the simulation.Ligand RMSD indicated fluctuations in three ZINC  This corresponds with the MM-GBSA results (Table 1) that were used to estimate the binding free energy of the compounds where the binding interaction between the protein and ligand complexes is specified by the free energy binding.The crystal ligand was used as a control where its score was −90.2 kcal/mol.The top compounds picked had significantly higher binding energy to the DOR, with the lowest of the scores being −134.9 ± 11.3 kcal/mol.
Protein−Ligand Interactions of the Top Eight ZINC Compounds to the DOR.The residues involved in the top eight ZINC compounds binding to the protein receptor were analyzed as described in the Methods section with the Desmond SID.All interacting residues to the top eight ZINC compounds and crystal ligand were tabulated (Table 4).2D ligand Water bridges occurred in all eight ZINC compounds except for ZINC000078648574, which is the same compound that did not show hydrogen bonding to residue ASP128.Most of the compounds showed higher hydrophobic interactions and hydrogen bonding, in comparison to the crystal structure, leading to higher GI absorption.
Secondary Structure Maintained for the Top Eight ZINC Compound−Protein Complexes.As was done for the crystal DOR system, protein SSE plots were generated for the top 69 ZINC compound−protein complexes (Figures S17− S20).Only minor changes in alpha-helical content were observed in each of these complexes, indicating that the protein secondary structure was maintained.
RMSF Shows Fluctuation in Regions of the Protein with Respect to the Ligand.Protein RMSF was calculated for the top eight ZINC compounds (Figure S21).Higher values depicted by peaks are areas of the protein that fluctuate the most, such as the N-and C-terminals as well as the intra-and extracellular loops.ZINC000037556415 showed the greatest fluctuation at nearly 5 Å around residue 200, which is located in extracellular loop 2. All compounds showed similar fluctuation at the same residue positions between 1 and 3 Å.A small fl u c t u a t i o n f r o m Z I N C 0 0 0 0 7 8 6 4 8 5 7 4 a n d ZINC000020559278 was observed around residue 250.Using the crystal structure as a positive control, each complex showed the same or higher residual fluctuation throughout the simulations.
Network Analysis Revealed Communication among Different Regions of the DOR.Unweighted and weighted network models of the DPI-287/DOR system and the top eight ZINC compound systems were calculated, as described in the Methods section, to decipher the allosteric signal transmission pathway.The comparison of the unweighted network models between the systems shows that their connections are in good agreement with each other.Quantifying the correlation between the nodes in the weighted network model reveals similarities observed between the systems.All systems appear to have higher correlations between edges TM5 and TM6.Community network models were generated using the weighted network model, which grouped residues together that interact more frequently and more strongly than to residues in other communities.The systems that seem to communicate more similarly to the crystal complex are ZINC000020559278, Z I N C 0 0 0 0 7 8 5 1 5 8 6 4 , Z I N C 0 0 0 8 2 7 3 6 0 7 9 4 , Z I N C 0 0 0 0 7 8 6 4 8 5 7 4 , Z I N C 0 0 0 0 5 7 9 9 9 6 5 3 , a n d ZINC000006664413.The basis for the similarity is that the intracellular portions of TM5 and TM6 are in the same community represented by one single color.Another observation in most of the similar systems is that the extracellular portions of TM6 and TM7 belong to the same community while the same region in TM5 belongs to a different community.This trend occurred in ZINC000078515864, ZINC000827360794, and ZINC000078648574.The systems that had the most similarities in critical nodes with the crystal system are ZINC000020559278 with 11 of the same critical nodes, ZINC000078515864 with 13, and ZINC000078648574 with 11 (Table S4).When calculating the optimal paths of the Transmission and Toggle Switch for each of the top eight ZINC compounds, multiple systems showed great similarities to the crystal complex (Figure 11 and Table S5).ZINC000025329384 shared every residue as the crystal ligand with the Transmission Switch pathway and two out of the four residues of the Toggle Switch pathway.ZINC000020559278, ZINC000078515864, and ZINC000037556415 also had many similar residues to the crystal complex with the Transmission Switch and Toggle Switch.ZINC000827360794 had the most residues similar to the crystal with the Toggle Switch.

Selectivity Analysis on DOR, KOR, and MOR.
To determine the selectivity of the top eight ZINC compounds for the DOR, the Glide XP docking scores of the top eight ZINC compounds on the crystal conformation (CC) of DOR, KOR, and MOR were obtained using the docking protocol in our original VSW (Table 5, columns 4−6).The original docking scores (Table 5, column 3) on DOR with the specific conformation (CC, C1/the first populated MD conformation, C2/the second populated conformation) were compared these docking scores on DOR, KOR, and MOR with CC.Two notable features are identified: (1) The docking score on DOR with specific conformation (C1 or C2) of the same compounds (3− 8) is on average better than the docking on DOR with CC by −1.8 kcal/mol, supporting our virtual screening strategy using multiple conformations rather than just crystal conformation to enrich the top eight ZINC compounds.(2) The docking score on DOR with specific conformation (CC, C1, and C2) is on

■ DISCUSSION
The opioid epidemic has brought to light the need for better opioid alternatives in public health all around the world.The continued rise of opioid addiction and overdose will not stop until there are better therapeutic agents available.Researchers and scientists have discovered that the DOR shows potential in not only pain management but also neurological and psychiatric disorders.Agonists targeting the DOR are strongly believed to not display addictive or dependence properties, such as MOR agonists, having the potential to help combat the addictive opioid crisis today.Although previous studies have been done on the DOR, none to our best knowledge have ever utilized multiple conformations from MD simulations for HTVS.MD simulations are able to probe deeper into interactions and dynamics that happen in a system that cannot be obtained from a crystal structure alone.Sampling the conformations from the active state DOR and using them for HTVS offers an opportunity to find better potential agonists to be of therapeutic use.Running long MD simulations on each of the top ZINC compounds identified from our VSW allows us to determine if the compounds bound to the DOR will remain stable and their protein−ligand interactions.Here, we present the first study using the ensemble-based HTVS approach to discover potential agonists to target the DOR.Our VSW combining molecular docking, MD simulation, bioinformatics tools, and drug similarity search revealed the top 69 hits from the ZINC15 database (32 from crystal conformation, 11 from cluster conformation 1, and 26 from cluster conformation 2).MM-GBSA and SwissADME prediction analyses reduced the top 69 ZINC compounds to eight compounds.The predicted drug ADME properties helped to specify if the compounds were highly (GI) absorbent as well as BBB permeant.These were the properties most valued due to opioids coming in an oral form and the receptors being located in areas of the brain.Without these specific properties, the compounds could be ruled out.
The protein−ligand interaction further validated the choice of the top eight ZINC compounds as potential DOR agonists.The large dataset and extended HTVS method portray the best interactions between ligands to form a complex with a molecular target.The adverse effects and related articles of the selected compounds were checked through CAS SciFinder and PubChem, where the compounds showed no adverse effects.Out of the top eight ZINC compounds, ZINC000057999653 is patented to be useful for altering the lifespan of eukaryotic organisms.These results further validate the top hits to be potential agonists.
The use of dynamic network models based on MD simulations data has shown to be efficient in extracting correlated motions, allosteric signals, and signal transduction networks within complex systems.The correlated motions are thought to be linked to their activity, which is normally difficult to accurately distinguish through visualization of the MD simulations alone.In addition, the communities that are generated with the dynamic network are highly correlated and The original docking scores of DOR were compared to the receptors based on their crystal conformations, respectively.The total average docking scores of DOR, KOR, and MOR for the top eight ZINC compounds are shown as well.The total average docking scores of each receptor were calculated based on the positive values.CC: the crystal conformation; C1: the first populated conformation from the MD simulation of CC; C2: the second populated conformation from the MD simulation of CC.See Table S2 for the initial compound reference number.
−73 In our study, the dynamic network analysis aided in identifying similar communication systems between the crystal DPI-287 and top eight ZINC compound systems.These similarities became more apparent when comparing the weighted networks that base their correlated motion in the simulation trajectories to the connections between nodes in the networks.The community models show that compound systems of ZINC000078515864, ZINC000827360794, and ZINC000078648574 communicate the most similarly to the DPI-287 system.This analysis highlights how the structural differences of the ligands can have similar or different dynamics to the DOR due to the grouping of communities based on residues that interact strongly and frequently with one another.
Based on the results, the potential binding and agonistic effects of the top eight ZINC compounds are indicated.Ensemble-based structure HTVS is a useful approach to find potential molecules that could target the binding pocket of the DOR.After examining 17 million ZINC15 compounds using structure-based HTVS methods, the most potential hits were further examined by MD simulations followed by MMGBSA binding free energy analysis.Further experiments are required to validate these compounds as potent DOR agonists.Nonetheless, this study has the potential to assist in the efforts to aid in the opioid epidemic.Experimental studies can be conducted on these compounds to help in the efforts of opioid addiction.

■ CONCLUSIONS
The lack of opioid alternatives with non-addictive properties has prompted researchers to look for effective candidates on all opioid receptors (mu, kappa, delta), with delta showing promising effects.Computational studies are a cost-effective method to identify a new target of existing drugs.Since the DOR is an attractive target for therapeutic effects without significant adverse effects, we have exploited the conformational flexibility of the DOR to search novel ZINC15 compounds, which may aid in opioid addiction.1 μs MD simulation of the active crystal conformation of the DOR was used to generate the structure ensemble.Using the clustering method, two major conformations of the DOR were identified.A total of three conformations (crystal conformation and two MD generated conformations) were used in our VSW of zinc compounds (17 million), leading to 69 compounds with top Glide XP docking scores and diverse structures.To further validate these compounds, 200 ns MD simulations were carried out to check the stability of the docked complexes, and the predicted drug ADME properties were examined.Eight stable systems were identified using a combination of dynamic properties (RMSD, RMSF) and MM-GBSA binding free energy calculation.Each of the eight top compounds exhibited better binding energy than the crystal ligand and contained attractive drug properties.Although this study suggests these top eight compounds may serve as good drug candidates for the DOR, further experimental studies and risk−benefit assessment are needed to evaluate the therapeutic values of the mentioned novel compounds.This study shows flexibility modeling the DOR; using MD simulations is a powerful tool in identifying novel compounds that could potentially show no adverse effects and aid in the opioid epidemic.
(Table S1) Molecular properties calculated using QikProp; (Table S2) detailed information regarding various properties from the Glide XP docking and MD simulations of the top 32 ZINC compounds in reference to the crystal structure (CC), the top 11 ZINC compounds in reference to the first conformation from the crystal MD simulation (C1), and the top 26 ZINC compounds in reference to the second conformation from the crystal MD simulation (C2); (Table S3) the predicted ADME properties for the top 32 best compounds in reference to the crystal conformation (PDB ID: 6PT3), the top 11 compounds in reference to the first conformation from the crystal MD simulation (C1), and the top 26 compounds in reference to the second conformation from the crystal MD simulation (C2), including the crystal reference compound from the SwissADME server.; (Table S4) critical nodes identified from the Network Analysis for the DPI-287/DOR system and the top eight compound systems; (Table S5) the optimal path of the transmission and tyrosine toggle switch generated from the Network Analysis for the DPI-287/DOR system crystal complex and the top eight ZINC compounds; (Figure S1

Figure 2 .
Figure 2. Virtual screening workflow to identify top compounds for DOR from the ZINC15 drug-like library.

Figure 3 .
Figure 3. Simulation interaction diagrams after MD simulation of the DOR crystal structure (PDB ID: 6PT3).(A) RMSD plot from MD simulation of 1000 ns.(B) Protein secondary structure elements (SSE).Orange represents alpha helices, and blue represents beta strands.(C) RMSF graph of protein of the crystal complex structure.(D) 2D ligand−protein interaction diagram from the MD trajectory.The residue displayed interacted with ligand for at least 30% of the simulation time.(E) Protein−ligand contacts during MD simulations.Interaction fraction greater than 1 is because of multiple contacts on one residue.

Figure 4 .
Figure 4. Superimposition of the active crystal DOR structure (PDB ID: 6PT3) (cyan) in complex with crystal agonist DPI-287 (purple) with the most abundant conformations from the MD simulations with agonist DPI-287.(A) DOR crystal conformation superimposed with the first abundant conformation (yellow, 75%) in complex with agonist DPI-287 (orange) including ligand-only view.(B) DOR crystal conformation superimposed with the second abundant conformation (pink, 24%) in complex with agonist DPI-287 (green) including ligand-only view.(C) DOR crystal conformation superimposed with both abundant clusters and all three ligand poses.

Figure 6 .
Figure 6.Top two low-frequency vibrational modes from the normal mode analysis based on the MD simulation of the crystal conformation of the active DOR in complex with agonist DPI-287.(A) Mode 1 and (B) mode 2. The top 10 modes can be found in Figure S1 of the supporting document.

Figure 7 .
Figure 7. Network analysis.(A) Structural communities generated from the weighted dynamic network model separated by color.(B) The critical nodes and weighted edges (red) shown for the crystal DPI-287 system with the node for the ligand represented in yellow (top).(C, D) Optimal signal transduction pathway of the Transmission switch (pink nodes) and Toggle switch (orange nodes) starting from the ligand (bottom).
compounds (ZINC000025329384, ZINC000037556415, and ZINC000078648574), whereas the remaining five compounds remained more stable throughout the simulation.MD Simulations Improved the Binding Pose of the Top Eight ZINC Compounds.Comparison between the Glide XP docking pose and the MD simulation pose was done for the top eight ZINC compounds (Figure 9).The same comparison was made for the remaining 61 ZINC compounds in the supporting document (Figures S10−S12).The simulation can significantly alter the ligand original bound conformation to optimize their interactions with the receptor.The simulation thus improved the binding pose of each of the top compounds.

Figure 8 .
Figure 8. Cα RMSD of the top eight ZINC compounds during 200 ns MD simulation in reference to the crystal active DOR conformation (PDB ID: 6PT3).

Figure 9 .
Figure 9.Comparison of Glide XP docking pose (blue) and MD simulation pose (red) for the top eight ZINC compounds.Receptor is in surface representation (gray).

Figure 10 .
Figure 10.2D ligand interaction diagrams for the top eight ZINC compounds with the crystal DOR receptor binding site.

Figure 11 .
Figure 11.Optimal signal transduction pathway of the Transmission Switch (pink nodes) and the Toggle Switch (orange nodes) starting from the ligand of each of the top eight ZINC compounds.
) top 10 low vibrational modes from the normal mode analysis (PCA) based on the active conformation DOR agonist DPI-287 system; (Figure S2) top 10 FDA-approved drug compounds based on their docking scores for the crystal conformation of the DOR; (Figure S3) top 32 ZINC compounds for the DOR crystal conformation including the crystal reference compound with ZINC ID, structure, SMILE code, docking score, number of STARs (indicator for "drug-likeness"), cluster IDs (ligand similarity clustering based on Canvas), and centroid (measures distances to the arithmetic means of clusters; (Figure S4) top 11 ZINC compounds for the first representative conformation of the DOR MD simulation structure conformation including the crystal reference compound with ZINC ID, structure, SMILE code, docking score, number of STARs (indicator for "druglikeness"), cluster IDs (ligand similarity clustering based on Canvas), and centroid (measures distances to the arithmetic means of clusters); (Figure S5) top 26 ZINC compounds for the second representative conformation of the DOR MD simulation structure conformation including the crystal reference compound with ZINC ID, structure, SMILE code, docking score, number of STARs (indicator for "drug-likeness"), cluster IDs (ligand similarity clustering based on Canvas), and centroid (measures distances to the arithmetic means of clusters); (Figure S6) comparison between the DOR crystal structure (PDB ID: 6PT3) (cyan) and the docked complex of the top 32 ZINC compounds (gray) in the side view and ligand view with the 2D chemical structure of the ZINC compounds; (Figure S7) comparison between the DOR crystal structure (PDB ID: 6PT3) (cyan) and the docked complex of the top 11 ZINC compounds (gray) of the first representative structure in the side view and ligand view with the 2D chemical structure of the ZINC compounds; (Figure S8) comparison between the DOR crystal structure (PDB ID: 6PT3) (cyan) and the docked complex of the top 26 ZINC compounds (gray) of the second representative structure in the side view and ligand view with the 2D chemical structure of the ZINC compounds; (Figure S9) Cα RMSD of the top ZINC compounds during 200 ns MD simulation in reference to the crystal active DOR conformation (PDB ID: 6PT3); (Figure S10) comparison of the top 32 ZINC compounds for the crystal conformation (PDB ID: 6PT3) in the docked pose (blue) and the MD simulation pose (red) with the DOR in surface representation (gray); (Figure S11) comparison of the top 11 ZINC compounds for the first representative structure from the MD crystal conformation (PDB ID: 6PT3) in the docked pose (blue) and the MD simulation pose (red) with the DOR in surface representation (gray); (Figure S12) comparison of the top 26 ZINC compounds for the second representative structure from the MD crystal conformation (PDB ID: 6PT3) in the docked pose (blue) and the MD simulation pose (red) with the DOR in surface representation (gray); (Figure S13) protein−ligand contacts during MD simulations for the top eight ZINC compounds; (Figure S14) protein− ligand contacts during MD simulations for top 32 ZINC compounds for the crystal conformation (PDB ID: 6PT3); (Figure S15) protein−ligand contacts during MD simulations for top 11 ZINC compounds for the first representative structure from the MD of the crystal conformation (PDB ID: 6PT3); (Figure S16) protein− ligand contacts during MD simulations for top 26 ZINC compounds for the second representative structure from the MD of the crystal conformation (PDB ID: 6PT3); (Figure S17) protein secondary structure elements (SSE) of the receptor in complex with the top eight ZINC compounds; (Figure S18) protein secondary structure elements for the top 32 ZINC compounds for the crystal conformation (PDB ID: 6PT3); (Figure S19) protein secondary structure elements for the top 11 ZINC compounds for the first representative structure from the crystal conformation (PDB ID: 6PT3); (Figure S20) protein secondary structure elements for the top 26 ZINC compounds for the second representative structure from the crystal conformation (PDB ID: 6PT3); (Figure S21) the Cα root mean square fluctuation (RMSF) of the receptor in complex with the top eight ZINC compounds; (Figure S22) the predicted ADME properties for the top 32 ZINC compounds based on the crystal conformation (PDB ID: 6PT3) including the reference compound from the SwissADME server (PDF) ■ AUTHOR INFORMATION Corresponding Author Chun Wu − Department of Molecular & Cellular Biosciences, College of Science and Mathematics, Rowan University, Glassboro, New Jersey 08028, United States; orcid.org/0000-0002-0176-3873; Email: wuc@rowan.edu

Table 1 .
Comparison between MOR, DOR, and KOR a a Note: MOR is mu-opioid receptor, DOR is delta-opioid receptor, and KOR is kappa-opioid receptor.

Table 2 .
Various Properties of the Top Eight ZINC Compounds Identified from our Virtual Screening Workflow a CC: the crystal conformation; C1: the first populated conformation from the MD simulation of CC; C2: the second populated conformation from the MD simulation of CC.See TableS2for the initial compound reference numbers.VDW is van der Waals interaction, and ELE is electrostatic interaction.b Based on the snapshots from the last 20 ns of simulation. a

Table 3 .
The Predicted Pharmacokinetics ADME Properties for Top Eight ZINC Compounds by SwissSimilarity Server a a Note: GI is gastrointestinal, and BBB is the blood−brain barrier.

Table 4 .
Protein-Ligand Interactions during MD Simulations for the Top Eight ZINC Compounds from the MD Simulations a average better than the docking score on KOR and MOR with CC by −1.3 and −2.4 kcal/mol, respectively.This clearly shows the good selectivity of our top eight ZINC compounds to DOR rather than KOR and MOR, if the binding entropy is comparable for these complexes.Only in two cases (ZINC000078515864 and ZINC000827360794) are the docking scores on DOR

Table 5 .
Glide XP Docking Scores of the Crystal Ligands DPI-287, MP1104, and BU27 and the Top Eight ZINC Compounds on the DOR, KOR, and MOR a