μ Opioid receptor: novel antagonists and structural modeling

The μ opioid receptor (MOR) is a prominent member of the G protein-coupled receptor family and the molecular target of morphine and other opioid drugs. Despite the long tradition of MOR-targeting drugs, still little is known about the ligand-receptor interactions and structure-function relationships underlying the distinct biological effects upon receptor activation or inhibition. With the resolved crystal structure of the β-funaltrexamine-MOR complex, we aimed at the discovery of novel agonists and antagonists using virtual screening tools, i.e. docking, pharmacophore- and shape-based modeling. We suggest important molecular interactions, which active molecules share and distinguish agonists and antagonists. These results allowed for the generation of theoretically validated in silico workflows that were employed for prospective virtual screening. Out of 18 virtual hits evaluated in in vitro pharmacological assays, three displayed antagonist activity and the most active compound significantly inhibited morphine-induced antinociception. The new identified chemotypes hold promise for further development into neurochemical tools for studying the MOR or as potential therapeutic lead candidates.

Supplementary Section S4. Evaluation and discussion of the applied virtual screening methods.
As results of different computational approaches, two active ligands at the MOR, 3 and 4, were identified with pharmacophore modeling one molecule, 5, was identified with shape-based modeling, and no active compound was among the top-ranked docking hits. It is well-known that docking is often limited by scoring and the correct ranking of compounds. 21 Since we selected our test compounds exclusively based on the GoldScore and the predicted interaction patterns, and did not include any manual inspection to avoid bias, this outcome was not surprising. In our recent study on cyclooxygenases (COX), 22 also docking performed modestly compared to the other methods when it came to enrichment of active compounds at the top of the hitlist. However, we identified one active molecule with a totally distinct chemical structure from all so far known COX-inhibitors. 22 This compound was ranked within the top-20 at position 16, so since we could only test the top-two and top-four ranked compounds, respectively, we might have missed actually active compounds ranked lower. For projects that aim at identifying structurally novel compounds, and therefore, a higher number of false positive hits can be accepted, docking might still offer advantages.
Although docking did not directly contribute to the identification of novel MOR ligands, still, the docking-based investigations represented the crucial pre-requisite for the generation of the pharmacophore-and shape-based models. The detailed analysis of the docking poses of agonists, antagonists, and inactive compounds, allowed for further insights into the mode of action, and without that knowledge, the previous model generation had failed. The agonist pharmacophore model generation performed very well, especially when a lot of XVols were added. However, compounds that are not allowed to protrude into a certain area of the MOR binding site also cannot form interactions in this site, which might explain the importance of the XVols for the agonist modeling. Despite the good performance of the models in the theoretical validation, they did not S17 retrieve any active molecules in the prospective part of the study. Again, this might be attributed to the low amount of hits selected for biological testing.
The MOR antagonist pharmacophore modeling appeared to be much more complex, as two of the models could not discriminate between agonists and antagonists in the theoretical validation. These two models may therefore rather represent general MOR ligand models than specific antagonist models.
Those two models also predicted the two top-ranked hits selected for biological testing. In detail, 3 mapped model pm-ant-lig-model-3, while 4 matched model pm-ant-lig-model-4. Interestingly, both molecules were found to be active at the MOR in the biological testing. The performance of the antagonist models is especially notable, considering the limited amount of experimental data.
Both agonist and antagonist shape-based modeling strategies proved to be very effective in the theoretical validation, and in addition the antagonist shape-models could discriminate between agonists and antagonists. This can be explained by the different definition of a "hit" applied in ROCS 23,24 and LigandScout. 25 In LigandScout, a compound needs to map all feature included in the model to be considered as hit. Subsequently, a pharmacophore fit value is assigned, that determines how well the compound maps the model. In ROCS, a score is assigned for every compound of the screening database, and the number of the final output results is user-defined. However, that implies that a compound does not have to match all features included in the model to be ranked. According to our hypothesis, there might be multiple options for antagonists to form additional polar interactions and prevent the activation of the receptor. It might be difficult to cover all these potential interactions, and the combinations thereof, in a few pharmacophore models as every compound needs to fulfill the whole hypothesis to be considered as hit. However, every compound that matches an additional polar color feature in a ROCS model (potential antagonist) gets ranked at a higher position as compounds that do not (potential agonists).
Therefore, ROCS models may have advantages especially for projects, where it is not clear which polar interactions might actually be crucial for binding of a ligand and which ones can be neglected, or for S18 projects, where many very similar binding modes can occur. In addition, also the activity of compound 5 predicted with the ROCS agonist model shape-ag-model-2 could be confirmed in the biological assessments. However, this compound displayed antagonistic activity.
To further examine the hit lists derived from the different applied methods we analyzed whether topranked compounds from one method were also predicted by other methods above the activity cut-off.
Several of these compounds were consensus hits, however, none of them was active in the biological evaluation. Specifically, T13 (predicted as agonist by docking) also mapped an antagonist pharmacophore model, the ROCS agonist hit T6 also mapped an agonist pharmacophore model and was predicted as agonist by docking with a GoldScore of 52.8, the agonist pharmacophore modeling hit T4 also fulfilled the docking criteria for agonists and retrieved a GoldScore of 44.4, and the agonist pharmacophore hit T2 was predicted as antagonist in docking since it formed additional interactions with Gln124 at a GoldScore of 57.9. This finding is in contrast to our previous study on COX, 22 where we could observe a strong correlation relationship between the number of methods that predicted a compound to be active and the actual activity in the experimental setting. Nevertheless, also in the same study 22 we have observed some cases where methods contributed to the identification of novel active molecules alone, and such findings were also reported by others. 26 The structural diversity compared to the known active ligands at the MOR may also explain the limited performance observed for the applied bioactivity profiling tools. The investigation of different ChEMBL 27 datasets in the course of the CYP study 28 highlighted the dependence of the 2D similaritybased tool SEA 17 on the underlying dataset of known active compounds. 28 An increasing number of active compounds was correctly predicted with increasing version of the ChEMBL database. Also, the experimental data available for the prediction increased in higher ChEMBL versions, thereby suggesting that the performance of the tool can be improved with the addition of further biological activity data.

S19
Since the ligands identified in the present work represent novel scaffolds, no active compounds might be included in the dataset that cover also these structural classes.
In the case of PharmMapper, 19 where also no molecule was mapped to a MOR pharmacophore model, we assume that this target may not be included in the pharmacophore model collection since the crystal structure was published not until recently. 1 To further evaluate this aspect, we profiled the co-crystallized MOR ligand β-FNA, which should under any circumstances map the model. However, also β-FNA was not predicted to interact with the MOR. Since the underlying reason remains speculative, we highly encourage the distributors of bioactivity profiling tools to provide information about the targets that are represented by their software. This novel interaction disrupts the cation-pi interaction between K303 and W318 in the 4DKL structure, leading to a shift of W318 and H319 in 5C1M. AS a consequence, the hydrogen bond between H319 and Y128 in 4DKL is disrupted as well and Y75 adapts its conformation to form a hydrogen bond with H319.

S29
14-O-benzylnaltrexone 6 14-O-benzylnaltrexone derivative 6 14-O-benzylnaltrexone derivative 6 14-O-benzylnaltrexone derivative 6 Table S4. Smiles codes of the molecules in the "inactives" dataset.   Pharmacophore modeling deploys the spatial arrangement of electrochemical features to represent the binding mode of ligands for a specific target 46 and to identify novel molecules that display similar feature patterns. Shape-based screening, as the name already indicates, uses the 3D shape of a known active compound for model generation and filters compounds with similar shapes out of large compound databases. The shape-models can further be modified via addition of chemical features, which define again electrochemical properties. 47 2D similarity-based methods compare the 2D structures of a set of known active compounds to a query molecule and employ statistical methods, based on the structural similarity, to calculate the probability that the compound of interest may also interact with the target.

Pharmacophore modeling
All pharmacophore modeling studies were performed with LigandScout version 3.1 25 . This program allows for both structure-and ligand-based modeling. In case a structure-based approach is employed, the program automatically calculates the interaction pattern present in the crystal complex. This initial model can then be refined with the "actives" and "inactives" dataset. In the course of ligand-based modeling, conformations of two or more known active are aligned, and common pharmacophore features are extracted. 25 Again, this initial model can then be manually refined. As default, LigandScout uses HBD, HBA, metal binding (M), PI, negatively (NI) ionizable, Aro, and H features. In addition, steric S37 constraints (XVols) can be added to a model to prevent mapping of compounds that would probably clash with the binding site. 25 All datasets and the Maybridge database were converted to 3D multiconformational databases using the Generate Library protocol (OMEGA version 2.3.3 [48][49][50] ) implemented in LigandScout. A maximum number of 50 conformers were generated for every molecule. Only compounds mapping all features included in the model were considered as hits.
After the prospective virtual screening, an additional molecular weight filter was applied. The molecular weight of the agonists or antagonists was calculated with Discovery Studio, and the average ± the standard deviation was used as filter. For agonists, the molecular weight ranged from 316 -473, whereas the antagonist hitlist was filtered for compounds with a molecular weight from 344 -506.
Finally, all hits were ranked according to the relative pharmacophore fit value, and the top 4-ranked agonists, and the top-2 ranked antagonists were selected for biological testing.

Shape-based modeling
The software program vROCS version 3.0.0 23,24 was employed for shape-based modeling. This program applies atom-centered Gaussian functions to calculate the shape-overlap of a query molecule and the model. In addition, so-called color features can be added. As default, HBD, HBA, R, H, anionic (A), and C colour features are available. Automatically generated models derived from the co-crystallized ligand or one low-energy conformation of the query molecules calculated with OMEGA 2.3.2 49-51 were further optimized with the "actives" and "inactives" dataset. The default Implicit-Mills-Dean force-field was applied for virtual screening, and the ComboScore was selected as scoring function.
The hit lists were filtered with the same molecular weight filter applied in the pharmacophore-based studies. In the course of a ROCS-based virtual screening run, a score is assigned for all compounds up to a user-defined number, no matter if the database compounds match all features or not. To discriminate S38 between active and inactive compounds, an activity cut-off was defined for every model during the theoretical model validation. These values can be quite distinct, as the probability that compounds match the majority of the features largely depends on the number and composition of the features. A model with a low number of features therefore more often retrieves compounds with a very high ComboScore compared to a complex model with multiple color features. To account for this bias, we normalized the retrieved ComboScore with the validation cut-off, and used this normalized ComboScore for ranking of the compounds in the prospective screening. Finally, the four top-ranked compounds from the agonist hitlist and the top-2 ranked compounds from the antagonist hitlist were selected for biological testing.
All datasets and the Maybridge database were converted to an oeb.gz file using OMEGA 2.3.2 and the default settings, leading to a maximum number of 200 conformers per molecules.

Docking
All docking studies were performed with GOLD version 5.2. 52 This program uses a genetic algorithm to calculate up to ten docking poses per input-ligand. The GoldScore was selected as fitness score, which takes hydrogen bonding, ligand internal strains, and steric aspects of the receptor-ligand complex into account. 53 The crystal structure of the β-FNA -MOR complex (PDB-entry 4DKL 1 ) was used for the docking studies. In the course of the protein preparation, hydrogens were added, and all water molecules except for 718 and 719 were deleted. The remaining two water molecules were set to "toggle and spin". This allowed the program to automatically decide whether or not a water molecule is included during docking, and to optimize the orientation of the water molecule. The area of 6 Å around the co-crystallized ligand was defined as binding site, and a protein hydrogen bond with the phenolic hydroxyl-group of Tyr148 was defined as constraint. Using these parameters, β-FNA was re-docked into the empty binding pocket of the MOR with an RMSD-value of 1.4330 Å. This RMSD-value might appear high at first sight, S39 however, a covalently attached ligand was docked using a "non-covalent" docking protocol. We decided for a non-covalent docking procedure, because the majority of known ligands and all of the compounds in our datasets were competitive binders. The double bond involved in the chemical reaction is about 3.6 Å apart from Lys233 in the docking pose, suggesting an appropriate position for the reaction. One lowenergy conformation per molecule calculated with OMEGA 2.3.2 49-51 served as input file for docking.
Since antagonists tended to have a higher GoldScore, two different cut-off values were defined for agonists and antagonists. Every agonist with a GoldScore ≥ 44.0 and antagonists with a GoldScore of ≥ 48.0 were considered as active. The interaction patterns of the highest ranked docking pose served as discrimination between agonists and antagonists.
In the prospective screening, a maximum number of ten docking poses for each of the ~ 52,000 compounds in the Maybridge database was derived. To limit the amount of data, the results were filtered according to the GoldScore, where all results above a GoldScore of 65.0 were kept. In addition this file was filtered one time with the agonist-molecular weight filter for potential novel agonist, and one time with the antagonist-molecular weight filter to identify novel antagonists, respectively. Both hitlists were then ranked according to the GoldScore, and the interaction patterns of the docked compounds were visually inspected. The four and two compounds with the highest GoldScore that fulfilled the required interaction patterns were selected for the biological evaluation as MOR agonists and antagonists, respectively.

Pharmacophore-based profiling
PharmMapper. The PharmMapper model collection contains more than 7,000 structure-based pharmacophore models covering approximately 1,600 drug targets. Although the models have been automatically generated with LigandScout, an in-house virtual screening algorithm was employed. 19 One single sd-file was uploaded for every compound, and a maximum number of 300 conformers was generated. The number of results was limited to the 300 best-matching pharmacophore models. To ensure a sufficient restrictivity, only models with ≥ 6 pharmacophore features were considered as hits.

S41
PharmaDB. In our previous studies, 22,28 we also included the commercial pharmacophore model collection PharmaDB, which was implemented in Discovery Studio. The models were automatically generated from target-ligand complexes deposited in the Protein Data Bank, and the restrictivity of the created models was estimated with a trained genetic function approximation model. A number of 3 to 6 pharmacophore features was allowed for every model, and a maximum number or ten models per complex was selected for the collection. 20 Unfortunately, this model collection does not contain any MOR models, however, for the sake of completeness we considered it as an added-value to the present study.
Hardware specification. All processes and predictions were performed on a multi-core workstation with 2.4+ GHz, 8 GB of RAM, a 1+TB fast mass storage, and a NVIDIA graphical processing unit. All programs run on the Windows 7 platform.