Simulating Protein Structure to Support Clinical Decisions: Predicting the Ligand Specic Impact of Genomic Variation on Opioid Binding

Background: Predicting the impact of missense protein variants on drug binding would have a widespread implications on the practice of genomic medicine, including matching a molecular therapy and dosage to an individual’s genome sequence. Genetic variation is widespread within G-protein-coupled receptors, which can affect overall structure and conformation of the receptors. These structural changes in turn impact ligand binding interactions, which may change the overall dosage requirements for target drugs. In this study, we used molecular docking simulations to explore the effect of missense variants on opioid drug binding anity to the opioid receptor mu 1 (OPMR1). Methods: Using high-throughput, in silico docking simulations, the binding interactions of 27 opioid drugs to naturally occurring variants in opioid receptor mu 1 (OPRM1) were used to predict changes to ligand binding anity. The binding energy of the small molecules to the wild-type receptor was compared to an experimentally derived inhibitory constant (Ki) for validation, and the variant-induced disruptions in variant:drug interactions used to predict the impact on the effective drug dosage. Results: The binding energies for each drug-variant receptor pair relative to the wildtype receptor and drug showed trends across drugs, with some variants showing enhancing (238I, 302I) or diminishing (235M, 235N) effects on binding anity. The 153V variant showed increased binding anity for certain drugs, and decreased anity for others. The simulation results correlated well with experimentally derived inhibitory constants (R 2 = 0.69), and an exponential regression model revealed how changes in relative binding energy between wildtype and variant structures predict changes to Ki. Conclusions: The simulation results illustrate the potential for integrating genetic variation into the process of development of small molecule therapies to support genomic-driven medicine. Depending on the drug and location, amino acid variation can either increase or decrease the strength of the molecular interactions and should be considered when determining drug dosage. The of variation the cost of experimental

inhibitory constants (R 2 = 0.69), and an exponential regression model revealed how changes in relative binding energy between wildtype and variant structures predict changes to Ki.
Conclusions: The simulation results illustrate the potential for integrating genetic variation into the process of development of small molecule therapies to support genomic-driven medicine. Depending on the drug and location, amino acid variation can either increase or decrease the strength of the molecular interactions and should be considered when determining drug dosage. The scale of variation and the cost of experimental characterization underscores the potential for accurate simulation based methods to inform clinical decisions.

Background
Genomic pro ling studies of healthy individuals have revealed an enormous amount of naturally occurring variation, and understanding the functional signi cance of these differences presents a major challenge to genomic medicine. This is especially true for missense variants, which presumably persist in the population because the underlying change to the protein structure only minimally impacts the wildtype function. However, when these variants occur in the functional domains of small molecule drug targets it can alter the drug response (1).
Predicting the impact of genomic variation on ligand binding would help inform clinical decisions when selecting an appropriate molecular intervention and dosage. This is particularly important when adverse drug responses can lead to a public health crisis, as with opioid misuse (2). While the driving factors underlying tolerance, addiction, and abuse are complex, the responsible administration of these powerful interventions would bene t from insight into how speci c drug:target interactions will be altered by missense mutations. Previous studies have shown that overprescription of opioids for analgesia can lead to overdependence and addiction (3,4), thus a better understanding of how genetic variation affects receptor binding sensitivity can lead to more appropriate dosage during treatment.
The genetic variation in the mu-opioid receptor gene (OPMR1) has been well-characterized, and many studies have shown an association between certain variants and opioid dependence (5-7). For certain highly-prevalent variants, the functional impact of the mutation on receptor binding a nity has been studied (8,9). However, due to the logistical constraints involved in in vivo functional assays, many minor variants located in the receptor binding pocket have not yet been studied.
The combinatorial explosion that results from the limitless pool of potential variants and continuously expanding list of drug candidates necessitates the use of a computational method to inform clinical decision making. Molecular docking is a common approach to modelling drug:target interactions, ranking various binding modes and predicting the binding a nity of a protein-ligand pair. Docking simulations have revealed information about receptor binding site conformations (10), and recently, identi ed novel drugs for the mu-opioid receptor (11)(12)(13). However, molecular docking studies focus primarily on the wild-type receptor, failing to capture the effect of functional variants on receptor binding a nity.
Our previous work highlighted the translational potential of a simulation-based approach for predicting the impact of drug-resistant somatic variation on available kinase inhibitors (14), leading to an improved work ow (SNP2SIM (15)). Using the SNP2SIM platform, we show the potential of a high-throughput in silico functional assay in studying G-protein-coupled receptors (GPCRs), in particular the opioid receptor mu 1 (OPRM1).

Methods
The initial structure of active OPRM1(homology model based on PDB: 5C1M) was downloaded from the G protein-coupled receptor database (GPCRdb.org). Of the dozens of natural variants annotated in the database (16), the six variants that mapped to the binding interface were used to generate mutant protein structures using DeepView (17) (Fig. 1). The opioid drug library was built from 27 opioid structures downloaded from DrugBank (18).
Binding simulations were run using the drugSearch module of the SNP2SIM (15), a Python-based work ow built to enable large scale computational simulations of variant protein structures. The drugSearch module utilizes AutoDock Vina (19) to perform docking simulations, automatically preparing the necessary input PDBQT les (using AutoDockTools) for both ligand and protein variant, and allowing for high-throughput docking experiments. For each variant-drug pair (189 combinations total), 10 independent simulations were performed using the default parameters and exhaustiveness of 50, and the binding energy of the highest ranked binding mode was averaged across trials to determine the predicted binding energy. For each simulation, the amino acid residues in the binding pocket were considered exible during docking, and the overall search space was determined using the crystal structure coordinates of these residues.
The varAnalysis module of SNP2SIM was used to calculate the change in binding energy (relative to the wild-type receptor-ligand pair), and visualize the results of the docking simulations. For a given drug, subtracting the wildtype binding energy from that of the variant allows for variant induced changes to drug interactions to be quanti ed. Positive values correspond to an increase in the interaction energy of the variant, which translates to a higher Ki for that particular drug.Using the relative binding energies for each variant also allows for a standardized comparison of the effect of each variant on the protein-ligand binding interaction. The varAnalysis module compiles the binding energy data from the outputs of the docking simulations, and calculates the mean relative binding energy across independent trials, as well as the variance. Each trial differs by the random seed used in Autodock Vina, and thus we can measure the range and uncertainty in binding energies for a variant protein-ligand pair. The varAnalysis module also outputs an R Shiny app for interactive visualization of the docking simulation results.
To validate the binding predictions to the wild-type receptor, the predicted binding a nities for each drugwildtype protein pair were compared to experimentally determined inhibitory constants (Ki) (20). Exponential regression was performed to model how in silico predicted binding energies relate to the experimental data (R 2 = 0.69). The regression model was generated without pentazocine, which was an outlier in Ki relative to drugs with similar binding results. This regression model was used to predict the change in Ki for each drug to a variant receptor, contextualizing the magnitude of change in binding energy in terms of a biological constant.

Results
Comparing the predicted binding a nity of a variant to that of the wildtype provides a relative measure of the variant impact on drug binding (Fig. 2). Since lower energies correspond to greater thermodynamic stability, increases in the relative binding energy correlate to decreased binding a nity of the drug compared to the wildtype receptor simulations. By calculating the relative binding energy and normalizing by the binding a nity to the wild-type receptor, we can more fairly compare ligand-protein binding energy across ligands. The full simulation results can be found in Table S1, along with code to execute their visualization using an interactive R Shiny app.
In general, variants show a similar impact across the opioid library, enhancing (238I, 302I) or diminishing (235M, 235N) the interactions of OPRM1 with the ligand. However, depending on the ligand, the 153V variant predicted both diminished (methadone, fentanyl) and enhanced binding (oxycodone, oliceridine) compared to the wildtype structure. The 153V variant in particular showed larger magnitude shifts (with a smaller range) in binding energy for multiple ligands, compared to the relatively large variance in results for ligands and other variants.
The natural neurotransmitter endomorphin-1 did not show drastic change in binding energy across variant receptors, which seems reasonable given these variants are not associated with any adverse conditions and found in healthy individuals. Several other ligands display a similar pattern of robustness ( Fig. 3) across variant receptors. These ligands serve as possible candidates for widely prescribable pain treatment due to their robustness to patient genotype. Figure 4 shows the binding pro les for the naturally occurring endomorphin-1 and morphine, and the widely used synthetic drugs fentanyl and oxycodone. Endomorphin 1 exhibits a higher variance in binding energy across trials compared to the remaining 3 drugs. Fentanyl, morphine, and oxycodone all showed relatively decreased binding energy (higher binding a nity) in variants 238I and 302I, suggesting that lower dosages may achieve the same analgesic effect (while minimizing risk of overdependence). The 153V variant varies signi cantly among these drugs, suggesting that a more curated approach to drug prescription may be needed for patients with this particular missense variant.
The drug naloxone is used widely to treat opioid overdose, and works by acting as a competitive antagonist for the opioid receptor. The drug is also often used in combination with other opioids to minimize the risk of overdependence, though the correct dosage in these situations still contains uncertainty (21). The docking simulation results indicate that naloxone may be robust to variants in the mu opioid receptor, though the relative binding energy decreased slightly in the 302I variant. However, no variant showed an increase in relative binding energy, indicating that individualized dosage might not be necessary when administering naloxone to prevent an overdose, which is often done in emergency situations where individual genotype information is not available.
The docking simulation results for each drug to the wild-type receptor were compared to experimentally determined Ki inhibitory constants (20), and an exponential regression model (Fig. 5) was built to predict the magnitude of change in Ki given the predicted change in binding energy for a drug-variant pair, relative to the wildtype receptor binding energy. Figure 6 shows the predicted range of Ki values for each variant-drug combination, binned into three ranges: ≤ 1 nM, 1-100 nM, and > 100 nM. Similar patterns of enhancing and diminishing effects can be seen across variants, though most variants were robust to extreme changes in predicted Ki, relative to the wild-type. In certain drugs, multiple variants increase the Ki by an order of magnitude (hydromorphone, buprenorphine, oxymorphine, alfentanil), while in other drugs, variants have the opposite effect (codeine, propoxyphene, fentanyl). Of particular note, for the drug pentazocine, certain variants decreased the Ki by two orders of magnitude, indicating a highly signi cant change in drug binding a nity, that could affect the potency of the drug and require a lower overall dosage in patients with those speci c variants. Though the model excluded pentazocine as an outlier, the docking simulation results independently predict a trend of increasing the sensitivity for all variant receptors to the drug. The outlying Ki of pentazocine may also be explained by experimental variability and the fact that it is only a partial agonist for the mu opioid receptor (a characteristic which is not captured by simulated docking experiments) (20). The experimental error in determining inhibitory constants highlight the importance of more stringent validation experimental evidence for a clinically acceptable model.

Discussion
The emerging discipline of genomic medicine will rely on an individual's genome to guide treatment decisions. Knowledge of how missense variants affect drug speci c binding a nity can inform prioritization or dosage of a particular drug, tailoring drug therapy to the individual patient. This personalized medicine approach requires a better understanding of how variants affect drug binding interactions, which can be achieved through the development of robust, accurate models of variant protein-ligand docking.
These simulation results also provide an additional insight into structure-based drug design, speci cally in the eld of analgesics. For example, opioid drugs which show variant-binding pro les similar to the naturally occurring endomorphin-1 (hydrocodone, mitragynine, PZM21) may be good candidates to seed the next generation of targeted therapeutics, as they could be more widely used due to the robustness to mu opioid receptor variability. Furthermore, studying how the variant-binding pro les change with drug structure can inform opioid drug design, allowing for new synthetic opioid drugs catered to speci c receptor phenotypes.
This current study is limited to variants that interact directly with the bound ligand, despite many more OPRM1 variants curated in GPCRdb. Speci cally, these variants are present in the receptor region which was considered exible during docking simulations, and thus directly impact the binding a nity with a given ligand through changes to the molecular interactions. However, variants outside the binding pocket of the receptor can also affect ligand binding a nity by a different mechanism. Changes to the protein's primary structure will in uence conformational dynamics and propagate to alter the spatial organization of binding residues, affecting the positioning of the ligand and ligand-binding pocket interactions.
Therefore, a structure-based prediction would require accounting for these conformational dynamics. Considering the dynamics of the variant protein system will also improve the docking simulations for variants within the binding pocket, as more biologically relevant, variant-speci c structures can be considered instead of a single rigid structure derived from the open crystal structure of the wild-type protein. Additional modules of SNP2SIM enable these computationally-intensive simulations by running multiple MD simulations, and clustering the frames into representative "scaffold" structures. However, the current version of SNP2SIM only allows for high-throughput MD simulations and analysis for cytoplasmic or extracellular protein domains, and is not currently con gured to handle the extra requirements of membrane-embedded structures. The promising results of this preliminary study on variant receptor docking simulations warrant future study into how the conformational dynamics of the receptor affect opioid binding a nity, and extending this analysis to more missense variants in OPRM1.
Collectively, our results suggest that variation in the OPRM1 gene of the mOR should be considered when prescribing a particular therapeutic regimen, as changes to drug binding may have major implications for drug e cacy and potentially lead to misuse of opioids. Our simulation data shows that altered drug binding can translate to orders of magnitude differences in Ki, which can affect drug e cacy. Patients who are unable to achieve pain relief from opioids prescribed by their clinician may return to their prescriber for more opioids, continue to suffer in pain, or seek other opportunities for pain relief through illicit drugs. The implications of our ndings highlight the impact mOR variation can have on drug dosing regimens as some variants may increase the potential for opioid misuse, due to patient's inability to achieve adequate therapeutic relief from their prescriptions. These ndings lay the groundwork for a better understanding of how pain management can be controlled through pharmacogenomics.
The current guidelines for the clinical interpretation of variants (American College of Medical Genomics, Association for Molecular Pathology) justi ably do not allow for the sole use of computational evidence to make a strong assertion of a variant's functional impact, but this is largely in part to the uncertainty of methods that only provide a generalized or non-ligand speci c prediction. The ligand-dependent impact on binding observed in OPRM1 153V is lost when the prediction is categorical (pathogenic/disruptive or benign) or does not explicitly consider differences in the biochemical nature of individual ligands. SNP2SIM enables the protein-speci c characterization of variants of unknown signi cance using a highthroughput in silico functional assay, providing a more robust form of computational evidence to support their clinical interpretation and translation to clinical practice. The results of this study highlight this promising approach to studying missense variants and their effect on protein function and drug binding, and further underscore the potential for computational approaches to drug design, functional variant classi cation, and personalized medicine.

Conclusions
This study focused on the effects of multiple variants in the mu opioid receptor on drug binding using high-throughput docking simulations. The results of the simulations were validated using experimentally derived inhibitory constants, highlighting the potential for computational approaches to drug design and screening. Using a regression model, changes in relative binding energy can be contextualized in terms of nanomolar dosage, which could be used to inform clinical decisions of opioid dosage for analgesics, minimizing the risk of opioid over-dependence in patients. These results show promise towards the goal of identifying novel drugs that are effective for speci c variant proteins, and developing more individualized treatment options based on a patient's genotype. Abbreviations OPRM1 opioid receptor mu, MD:molecular dynamics Declarations Ethics approval and consent to participate Not applicable.