Comprehensive Signaling Profiles Reveal Unsuspected Functional Selectivity of δ-Opioid Receptor Agonists and Allow the Identification of Ligands with the Greatest Potential for Inducing Cyclase Superactivation.

Prolonged exposure to opioid receptor agonists triggers adaptations in the adenylyl cyclase (AC) pathway that lead to enhanced production of cyclic adenosine monophosphate (cAMP) upon withdrawal. This cellular phenomenon contributes to withdrawal symptoms, hyperalgesia and analgesic tolerance that interfere with clinical management of chronic pain syndromes. Since δ-opioid receptors (DOPrs) are a promising target for chronic pain management, we were interested in finding out if cell-based signaling profiles as generated for drug discovery purposes could inform us of the ligand potential to induce sensitization of the cyclase path. For this purpose, signaling of DOPr agonists was monitored at multiple effectors. The resulting signaling profiles revealed marked functional selectivity, particularly for Met-enkephalin (Met-ENK) whose signaling bias profile differed from those of synthetic ligands like SNC-80 and ARM390. Signaling diversity among ligands was systematized by clustering agonists according to similarities in E max and Log(τ) values for the different responses. The classification process revealed that the similarity in Gα/Gβγ, but not in β-arrestin (βarr), responses was correlated with the potential of Met-ENK, deltorphin II, (d-penicillamine2,5)-enkephalin (DPDPE), ARM390, and SNC-80 to enhance cAMP production, all of which required Ca2+ mobilization to produce this response. Moreover, superactivation by Met-ENK, which was the most-effective Ca2+ mobilizing agonist, required Gαi/o activation, availability of Gβγ subunits at the membrane, and activation of Ca2+ effectors such as calmodulin and protein kinase C (PKC). In contrast, superactivation by (N-(l-tyrosyl)-(3S)-1,2,3,4-tetrahydroisoquinoline-3-carbonyl)-l-phenylalanyl-l-phenylalanine (TIPP), which was set in a distinct category through clustering, required activation of Gαi/o subunits but was independent of the Gβγ dimer and Ca2+ mobilization, relying instead on Src and Raf-1 to induce this cellular adaptation.

The δ-opioid receptor (DOPr) is considered an attractive target for chronic pain management.[1] In effect, agonists that activate this receptor display analgesic efficacy in preclinical models of inflammatory,[2] neuropathic,[3−6] and cancer pain,[7] while their anxiolytic and antidepressant properties[8] provide a means of managing the distress associated with these chronic conditions.[9] Moreover, in comparison to μ-opioid receptor (MOPr) agonists, DOPr agonists have less potential for abuse[10,11] and physical dependence,[12,13] while displaying mitigated respiratory[10,14] and gastrointestinal side effects.[10,15] Nonetheless, development of acute[16,17] and chronic[18] analgesic tolerance remains a matter of concern although this side effect appears to develop in a ligand-specific manner.[6,19] Hence, a better characterization of the determinants of this specificity should help improve the development of more effective DOPr agonists for chronic pain management.

Among ligand-specific mechanisms of analgesic tolerance, the failure to support DOPr recycling has been identified as a reliable predictor of mitigated development of both acute[20] and chronic[6,19] tolerance to the analgesic actions of DOPr agonists. However, since analgesic tolerance is multifactorial,[21−24] mechanisms beyond receptor regulation should also be considered. Some of these adaptations such as sensitization of the cyclase pathway take place at the cellular level[25] and have been associated with the development of hyperalgesia[26] and analgesic tolerance[27−29] that develops upon repeated opioid administration. The use of transgenic murine models has provided insight into specific cyclase subtypes that contribute to the induction of analgesic tolerance by MOPr and DOPr agonists, especially the role of calcium-dependent adenylyl cyclases (ACI and ACIII)[29] and adenylyl cyclase V (ACV).[30] Moreover, the contribution of cyclase superactivation as a specific mechanism contributing to the development of analgesic tolerance has been documented.[27] In particular, repeated injection of morphine into the ventrolateral periaqueductal gray region induces sensitization of the cyclase path leading to enhanced neurotransmitter release by GABAergic interneurons in this structure and reduces analgesia by morphine.[27]

Signals that support cyclase superactivation and lead to enhanced cyclic adenosine monophosphate (cAMP) production following sustained exposure to DOPr agonists have been previously described.[31−33] On the other hand, the possibility that these adaptations may be triggered in a ligand-specific manner has not been addressed. The identification of drug candidates that have the potential to induce such unwanted adaptations is desirable, but doing so at early stages of the drug screening process remains a challenge, especially if ligands display functional selectivity.[34,35] We have previously shown that establishing similarities among signaling profiles at multiple signaling readouts allows one to group together GPCR ligands that will also share similar responses in more complex outcomes such as clinical side effects.[36] Hence, here, we wanted to determine whether classifying ligands according to signaling similarities at physiological DOPr effectors could be informative of their potential for inducing cyclase superactivation. To address this question, we first characterized the signaling profile of DOPr agonists at multiple downstream effectors and, then, we classified ligands according to similarities in logistic and operational parameters describing their signaling efficacy at these readouts. Finally, we assessed if these similarities were correlated with the ligands’ potential for sensitizing the cyclase pathway.[36]

Signaling profiles revealed unprecedented functional selectivity of DOPrs agonists, particularly for the endogenous ligand Met-enkephalin (Met-ENK) whose bias profile differed considerably from that of synthetic agonists such as SNC-80 and ARM390. Despite considerable functional selectivity of different ligands, the identification of their signaling similarities allowed us to successfully recognize (N-(l-tyrosyl)-(3S)-1,2,3,4-tetrahydroisoquinoline-3-carbonyl)-l-phenylalanyl-l-phenylalanine (TIPP) as a ligand with a unique response profile, while the estimates of the signaling similarity among Gα/Gβγ responses by the rest of ligands were directly correlated with their potential for inducing superactivation of the cAMP pathway via a Ca2+ sensitive mechanism.

Results and Discussion

The Endogenous Ligand Met-ENK Displays Unprecedented Signaling Diversity at DOPr Effectors

To characterize the signaling profiles of DOPr agonists, we used 12 different BRET-based biosensors that allowed us to monitor G-protein activation,[37] cAMP production,[38] Kir3 channel opening via Gβγ subunits,[39] and β-arrestin (βarr) recruitment in the absence or presence of different complements of GRKs.[36] Ca2+ mobilization was also monitored using obelin, a biosensor which produces luminescence upon Ca2+ binding.[40] Each response was monitored at the time of peak effect for each biosensor readout. Concentration response curves for effectors that directly interact with the receptor, such as G proteins and βarrs, are shown in Figure , and the curves for further downstream responses, such as cyclase inhibition, Ca2+ mobilization, and Kir3 channel activation, appear in Figure .

Figure 1

G protein activation (A) and βarr recruitment (B) by DOPr agonists were monitored using BRET-based biosensors (represented schematically). Results correspond to mean ± SEM, and the number of independent experiments for each readout and ligand are indicated in the figure. Responses elicited by different agonists were normalized to the maximal effect produced by Met-ENK, which was tested in all experimental runs. Curves were fit with the operational model and the logistic equation (curves shown correspond to logistic fits). Values for operational and logistic parameters are provided in Table S1 and are graphically summarized in Figure .

Figure 2

Ca2+ mobilization, adenylate cyclase inhibition, and Kir3 channel activation by DOPr agonists. Ca2+ mobilization (A), cAMP accumulation (B), and Kir3 channel activation (C) by DOPr agonists were monitored using BRET-based biosensors (represented schematically). Results correspond to mean ± SEM; the number of independent experiments for each readout and ligand are indicated in the figure. Responses elicited by different agonists were normalized to the maximal effect produced by Met-ENK, which was tested in all experimental runs. Curves were fit with the operational model and the logistic equation (curves shown correspond to logistic fits). Values for operational and logistic parameters are provided in Table S1 and are graphically summarized in Figure .

Figure 3

Graphic representation of operational and logistic parameters describing responses of DOPr agonists relative to DPDPE: Emax/EmaxDPDPE (A), ΔLog(τ) = Log(τ) – Log(τ)DPDPE (B), and transduction coefficients ΔLog(τ/KA) = Log(τ/KA) – Log(τ/KA) (C) were derived from concentration response curves in Figures and 2. Their values are given in Table S2, and here, they are shown as radial graphs.

The curves for G protein activation (Figure A) indicate that SNC-80 was a full agonist at these readouts and consistently produced greater G protein responses than Met-ENK. In turn, the maximal response of the endogenous ligand was marginally larger than that of the partial agonist (d-penicillamine2,5)-enkephalin (DPDPE) at Gαi1 and similar for both ligands at GαoA, Gαi2, and Gαz responses. In contrast, βarr recruitment curves for Met-ENK were practically superimposed with those produced with the full agonist SNC-80, while DPDPE remained a partial agonist across all readings (Figure B). These observations are consistent with Met-ENK being more effective at promoting βarr recruitment than G protein activation. The possibility that the endogenous ligand displayed biased responses at these readouts was verified by comparing transduction coefficients (details on curve fitting in Materials and Methods). The results of these comparisons using DPDPE as the standard ligand are presented in Table S2. Information provided in this table indicates that normalized transduction coefficients for βarr2 recruitment (±GRKs) by Met-ENK were ∼35–60-fold larger than those obtained for Gαi2 and Gαz stimulation and the transduction coefficient for βarr2 was 5-fold that of GαoA, whereas no significant difference between βarr and Gαi1 values was observed. Further analyses of Met-ENK effects across the remaining biosensors indicated that this agonist produced the largest mobilization of Ca2+ (Figure A) and was an effective inhibitor of cAMP production (Figure B). On the other hand, with the exception of TIPP that failed to induce Kir3 signaling, Met-ENK produced the smallest response at this readout (Figure C). Such diverse patterns of response across the different effectors resulted in additional bias in Met-ENK signaling, which included ∼10–80-fold preference in Ca2+ mobilization over activation of different G proteins, greater than 150-fold preference in favor of βarr signaling as compared to Kir3 channel response, and ∼250-fold more effective Ca2+ mobilization than Kir3 signaling (Table S2). Previous studies that examined biased responses by this endogenous opioid had reported a lack of bias at MORs[41] and had similarly identified that DOPr activation by Met-ENK favored βarr recruitment over Gα activation.[42] Here, we show substantial diversity of DOPr-mediated Met-ENK responses across a diversity of pathways.

Response Profiles of Other DOPr Agonists Differ from Those of the Endogenous Ligand Met-ENK

Concentration response curves for the nonendogenous DOPr agonists were also fit with logistic and operational equations to yield the corresponding parameters provided in Table S1 and graphically represented in Figure . Table S2 summarizes ΔΔLog(τ/KA) values (where KA is the affinity constant) for these same ligands, using DPDPE as the standard. From the information presented therein, it is evident that some of these nonendogenous agonists also displayed signaling versatility at different readouts but their bias profiles differed from those of Met-ENK. Indeed, if we consider Emax (Figure A) and Log(τ) (Figure B) values, we see that βarr recruitment and G protein/Kir3 activation by SNC-80 and ARM390 are consistently larger than those describing responses for DPDPE, while the same responses for TIPP are consistently smaller (or nonexistent). On the other hand, Met-ENK’s position relative to DPDPE is not the same for βarr recruitment vs G protein/Kir3 activation. Such differences result in bias for βarr vs G protein/Kir3 responses for Met-ENK but not for SNC-80 or its analogue ARM390[43] (Table S2). Met-ENK’s signaling bias in favor of βarr is also in contrast with the reported preferential activation of G protein by the novel DOPr agonist PN6047.[44]

Deltorphin II’s variations in Emax and Log(τ) values for βarr and G protein/Kir3 responses relative to DPDPE were less pronounced than for Met-ENK but more prominent than for SNC-80 and ARM90 (Figure A,B), pointing to some resemblance in βarr vs G protein signaling bias displayed by Met-ENK and deltorphin II (Table S2). Ca2+ signaling also differed among the two synthetic agonists and the endogenous ligand with Ca2+ mobilization by Met-ENK being the largest among all agonists tested. Such differences resulted in opposing signaling bias for Met-ENK as compared to SNC-80 and ARM390 in relation to Ca2+ mobilization vs G protein or Kir3 channel activation (Table S2). The marked potency of deltorphin II to inhibit cAMP production should also be noted. The overall difference in signaling profiles between naturally occurring peptides like Met-ENK and deltorphin II versus those of synthetic agonists is of interest and should be considered for the rational design of novel DOPr analgesics. In this sense, it is worth noting that in vivo administration of SNC-80 or ARM390 but not of deltorphin II produces analgesic tolerance upon repeated administration.[6,18] The distinct development of tolerance by deltorphin II and SNC-80 is partly determined by the unique ability of the peptide to support recycling.[6] In addition to their distinct postendocytic trafficking and different signaling profiles, synthetic agonists and peptidic ligands signal from different compartments. In particular, while peptides induce surface and endosomal signaling, synthetic agonists additionally initiate Golgi signaling.[45] Taken together, these distinctions could contribute to deltorphin II’s favorable tolerance profile in vivo.

Epileptogenic activity is another side effect that limits the clinical use of DOPr agonists. While SNC-80 produces seizures, ARM390 does not.[46] These ligands did not differ in the nature of signals produced but rather in their efficacy to produce them (Figure A), suggesting that highly efficacious ligands would be more prone to this side effect.

It is also intriguing that the endogenous ligand Met-ENK and synthetic agonists SNC-80/ARM390 displayed opposing preferences toward Ca2+ mobilization and Kir3 channel activation since for DOPr agonists both of these responses are mediated via activation of Gβγ effectors.[2,47] A possible explanation for these divergent actions is that Met-ENK induces Ca2+ mobilization via additional signals. Alternatively, these observations could also imply that the way Gβγ is released from the Gα “hotspot”[48,49] and/or the way in which the Gβγ dimer interacts with its downstream effectors PLCβ and Kir3[2,47] is ligand-specific. This latter possibility would require receptor and downstream effectors to be part of a complex that also includes the heterotrimeric G protein, a configuration that has been reported for Kir3 channels and GPCRs, including DOPr.[47]

DOPr Agonists Can Be Classified According to Signaling Similarities

The signaling profiles generated above revealed considerable signaling diversity of DOPr agonists at physiological effectors. This variety in signaling properties poses a challenge in terms of organizing information in the context of drug discovery, particularly to identify the signals that drive the pharmacological actions we seek to enhance or avoid in drug candidates.[34,35] Hence, once we had established that signaling diversity exists, we were interested in identifying what signals can be associated with the greatest potential for inducing cyclase superactivation, an adaptation known to interfere with opioid analgesia.[27−29] The question was addressed in two steps: First, we organized ligands by classifying them according to signaling similarities, and second, we assessed whether ligands within different signaling categories displayed a distinct potential for sensitizing forskolin-induced cAMP production.

Ligand classification was established using a previously described statistical method that uses operational and logistic parameters from multiple functional readouts to measure signaling similarities among ligands. We have previously shown that clustering ligands according to efficacy-related parameters (i.e., Log(τ), Emax) rather than additionally including affinity/potency information (Log/KA, EC50) allowed us to best associate signaling categories with responses of interest.[36] Thus, ligands were classified according to similarities in Emax and Log(τ) values (see Materials and Methods for details on curve fitting).

The measure of signaling similarity that is provided by the classification procedure corresponds to the frequency of coclustering of pairs of ligands across the iterative comparisons of their parameters built into the computational method.[36] These coclustering frequencies, which range from zero (a specific pair of ligands is never clustered together) to one (a specific pair of ligands is always grouped together), are then organized into a similarity matrix where each drug is assigned both a row and a column. Hence, if for example one wants to find out how similar SNC-80 and Met-ENK are according to the method applied, one identifies the column assigned to SNC-80 and the row assigned to Met-ENK and consults the frequency value in the corresponding cell of the similarity matrix. Figure A–C shows the similarity matrices for the indicated readouts in the form of a heatmap (where blue represents a frequency of one and yellow a frequency of zero coclustering). In these heatmaps, the rows and columns of the original similarity matrix were additionally rearranged by hierarchical clustering to highlight the groups of ligands with shared signaling properties as well as the corresponding clustering tree[36] (Figure A–C).

Figure 4

DOPr ligands can be clustered according to Log(τ) and Emax values, and estimates of similarity among these parameters are correlated with ligand potential to sensitize cAMP production. Heatmaps and dendrograms representing ligand similarity in Gβγ-mediated responses (Kir3/Ca2+) (A), Gα activation (B), and βarr recruitment (C). Yellow and blue, respectively, indicate ligands that never or always cluster together; coclustering frequencies resulting from the iterative comparison of parameters for the indicated pairs of ligands appear within each cell. The frequency of coclustering with Met-ENK was obtained for each ligand. These estimates of ligand similarity to the endogenous agonists for Gβγ-mediated signals (D), Gα activation (E), and βarr recruitment (F) were correlated to fosrkolin-induced cAMP production following exposure to Ca2+-mobilizing ligands. cAMP production was estimated from AUC values of forskolin concentration curves obtained in cells pre-exposed to the different DOPr agonists. r2 and p values are shown within the corresponding plots.

From the heatmaps in Figure , it is evident that signaling diversity among ligands remains manifest after clustering. For example, TIPP was set apart in a class by itself across the three types of readouts considered. The separation of this ligand from the rest is consistent with TIPP’s uniquely weak signaling efficacy, a characteristic that has been previously documented both in vivo(6) and in vitro.[2,19,50,51] Deltorphin II and DPDPE were similar when considering Gα- and Gβγ-mediated signals (Figure A,B) but less so with respect to βarr recruitment (Figure C). ARM390 and SNC-80 were in the same category in relation to Gβγ-mediated signals but not for the other readouts. In fact, SNC-80 was in a category by itself for Gα stimulation, representing the fact that it was the only ligand to consistently elicit maximal responses in all Gα subtypes assessed. In contrast, SNC-80 was clustered with Met-ENK for βarr recruitment since they both are similar in their capacity to induce maximal response in this readout. Met-ENK itself displayed varied positions across the different classifications, representing its unprecedented signaling diversity.

Similarities Among Gα/Gβγ Responses Are Correlated with Ligand Potential to Induce Cyclase Superactivation

Having established a signaling-based classification of the DOPr agonists, we next sought to determine whether ligands distributed into the different categories displayed distinct levels of cyclase superactivation. Hence, we started by evaluating the superactivation induced by the different DOPr agonists. To do so, HEK293 cells expressing the receptor and a biosensor monitoring cellular levels of cAMP were exposed to vehicle (DMSO, 0.1% (v/v)) or to the following agonists: Met-ENK, deltorphin II, DPDPE, TIPP, ARM390, or SNC-80 (10 μΜ; 8 h).[52] At the end of the treatment, the cells were washed and exposed to increasing concentrations of forskolin, a direct activator of cellular adenylyl cyclases[53] (Figure A,B). Preincubation with all of the agonists tested enhanced maximal cAMP production by forskolin. In addition, preincubation with Met-ENK, TIPP, and DPDPE also produced a significant increase in the potency (pEC50) of the forskolin response (Figure C,D). We did not identify any ligand that modified pEC50 values leaving Emax unchanged. The integration of the measures of cAMP levels across increasing concentrations of forskolin (area under the curve: AUC) provided an estimate of the overall increase in forskolin-driven second messenger production following exposure to the different agonists. Met-ENK produced the greatest sensitization to forskolin, and the rank order for the rest of the ligands was as follows: Met-ENK > TIPP ≅ DPDPE ≅ deltorphin II > ARM390 ≅ SNC-80 (Figure E). These differences in forskolin-induced cAMP production could not be simply attributed to residual DOPr activation after agonist washout. Indeed, basal cAMP levels after three washes were similar in cells pre-exposed to vehicle and cells exposed to different agonists (Figure S1), certifying no residual DOPr modulation of cAMP levels after washout.

Figure 5

Superactivation of the cyclase pathway by DOPr agonists. To find out how the cellular ability to produce cAMP was modified by sustained exposure to DOPr agonists, HEK293 cells expressing the receptor and a BRET biosensor that allows one to monitor cellular levels of cAMP were exposed to the indicated agonists (10 μM) for 8 h. At the end of treatment, cells were washed and exposed to increasing concentrations of forskolin (Fsk) before cAMP levels were monitored by BRET (A, B). Forskolin concentration response curves show cAMP accumulation in cells that were pre-exposed to vehicle (CTL) or to different agonists. Results represent mean ± SEM normalized to forskolin production in controls. The number of independent experiments for each condition are indicated in (C). Radial graphs representing logistic parameters (Emax, pEC50) and area under the curve (AUC) for forskolin responses following preincubation with DOPr agonists or vehicle (CTL). Dashed lines show 95% confidence interval limits. Emax is shown on a scale where CTL is 100%. AUC is shown in multiples of CTL AUC. Statistical comparisons between cAMP production observed in the CTL condition and following exposure to indicated agonists were done by verifying the lack of overlap for CI95 values (*), and CI99 values (**) (D). Rank ordering of the different agonists according to AUC for cAMP production is provided (drugs are deemed indistinguishable if the CI95 of their AUC overlap) (E).

Next, we sought to associate enhanced cAMP production by each ligand to the signaling categories that were generated with Gα-, Gβγ-, and βarr-derived parameters. To do so, measures of signaling similarity (i.e., frequency of coclustering) between each ligand and Met-ENK were retrieved from the corresponding similarity matrices/heatmaps (Figure A–C), and the retrieved frequency values were then correlated with AUC estimates of cAMP production (Figure D−F). Despite the noted diversity in signaling profiles, the measures of signaling similarity obtained by clustering were successfully correlated to cAMP levels. The strength of the association was most evident for the Gβγ-driven classification, where ligand similarity explained 98% of the variance associated with enhanced cAMP production by Ca2+-mobilizing ligands (Figure D) (r2 = 0.98; p = 0.001). Measures of Gα signaling similarity also explained a considerable proportion of the variance for these ligands (r2 = 0.89; p = 0.017; Figure E) as did coclustering of Gα and Gβγ parameters (r2 = 0.98; p = 0.001; Table S3). The classification strategy also recognized TIPP as a unique DOPr agonist. Consistent with TIPP being set apart from all other ligands, its inclusion in Gα and Gβγ correlations distorted the associations observed for Ca2+-mobilizing ligands (Table S3).

No association was observed between ligand similarity to induce in βarr responses and sensitization in cAMP production, independent of whether TIPP was included or not in the correlations (Figure E; Table S3). Furthermore, if instead of considering signaling similarities we considered parameters describing the curves for βarr recruitment or even AUC values derived from these curves, no association with cAMP levels was revealed either (Table S4), underscoring the independence of βarr signaling and sensitization of cAMP production for this set of ligands. No evident association was revealed between enhanced cAMP production and ligand-induced internalization of the receptor either (Figure S2). This lack of correlation is not unexpected given the independence between βarr recruitment and cAMP production following prolonged exposure to the different agonists.

It is well established that the time course of the different signals that results from the activation of a GPCR is quite distinct. Because of these distinct kinetics, the time frame of data acquisition may influence the estimation of bias magnitudes and overall assessment of functional selectivity when using bias as a descriptor.[54] Classification according to signaling profiles evaluates functional selectivity independent of bias measures. For the purpose of comparing and clustering together drugs with similar signaling profiles, the different signals were monitored at the time of peak response for each biosensor readout ensuring that maximal signaling efficacy at each biosensor was captured. The time to peak response takes place within minutes and does not represent the intricate evolution of signals over the 8 h time frame leading to cyclase adaptations. Nonetheless, efficacy measures taken at the time of the maximal response for each biosensor were successfully correlated with the protracted functional consequences of sustained exposure to an agonist. Similarly, when peak responses at MOPr were used to classify clinically available opioid analgesics, signaling similarities were highly predictive of the frequency of report for faecaloma,[36] a clinical manifestation that develops over days.

DOPr Agonists Require Gαi/o Activation to Sensitize the cAMP Pathway and They Do so via Gβγ-Ca2+-Calmdulin-Protein Kinase C or Src/Raf-1 Signaling Cascades

In the previous section, we observed that signaling similarities in Gα/Gβγ responses by Met-ENK, DPDPE, deltorphin II, ARM390, and SNC-80 were correlated to the potential of these ligands to sensitize cellular production of cAMP. We now want to evaluate whether the signals associated to superactivation by similarity are also mechanistically linked to this outcome. Of the similarity associations established above, Gβγ responses were the most highly correlated with enhanced cAMP production. Not only does Gβγ drive Ca2+ mobilization by DOPr agonists (see Gendron et al.[2] and references therein), but this second messenger modulates the activities of cyclases ACI and ACIII[55−59] that are highly expressed in HEK293 cells,[60] including those used in the present study (Figure S3). Moreover, Met-ENK, which was the ligand with the largest superactivation response (Figure A), was also the most effective at inducing Ca2+ mobilization (Figure A). Hence, we determined if mobilization of intracellular Ca2+ contributed to cAMP sensitization by the different agonists. To do so, the forskolin-promoted production of cAMP was compared in cells exposed or not to the intracellular Ca2+ chelator BAPTA (3 μM)[61] during agonist treatment (10 μM; 8 h). With the exception of TIPP, which failed to mobilize Ca2+ and was also set apart by the signaling classification, BAPTA interfered with the enhanced production of cAMP by all other ligands (Figure A). Also consistent with Gβγ signaling being a good predictor of cAMP sensitization, scavenging of the βγ dimer by overexpression of a membrane-bound version of the C-tail of GRK2[62] blocked cyclase sensitization by the prototypical Ca2+-mobilizing ligand Met-ENK but not by TIPP (Figure B). On the other hand, incubation with PTX (100 nM; 16 h) inhibited sensitization by both of these agonists (Figure C), as did the absence of the receptor or the introduction of the DOR antagonist naltrindole (10 μM) throughout the treatment phase (Figure S4), indicating that adaptations induced by Ca2+ mobilizing and nonmobilizing agonists both relied on the activation of the receptor and stimulation of downstream Gαi/o proteins. To further characterize the differential implication of Ca2+ in cyclase superactivation by Met-ENK and TIPP, we evaluated the contribution of calmodulin (CaM), a major effector in the modulation of Ca2+-sensitive cyclases.[56−59] Consistent with the fact that Ca2+ supported adaptations by Met-ENK but not TIPP (Figure A), the calmodulin blocker calmidazolium (10 μM) inhibited sensitization of cAMP production by the endogenous ligand (Figure A) but did not have an effect on the sensitization by TIPP (Figure E). Furthermore, chelerythrine (5 μM), which inhibits Ca2+-sensitive PKCs (protein kinase C),[63] further corroborated differences between the two ligands as it practically abolished superactivation by Met-ENK, leaving TIPP’s response unaffected (Figure B,F).

Figure 6

Ca2+ mobilization, Gβγ dimers, and Gαi/o signaling contribute to adenylate cyclase activation by DOPr agonists. HEK293 cells expressing DOPr and a BRET biosensor that allow one to monitor cellular levels of cAMP were exposed to the indicated agonists (10 μM; 8 h) in the presence (white) or absence (gray) of the Ca2+ chelator BAPTA (3 μM). At the end of the experiments, concentration response curves for forskolin were generated and AUC ± 95% CIs were calculated as in Figure . The number of independent experiments per condition are indicated at the bottom of each histogram bar. Statistical comparisons between responses observed in the presence or absence of BAPTA were done by verifying the absence of overlap for CI95 (*), CI99 (**), and CI99.9 (***) (A). The same cells as above were transiently transfected to express (white) or not (gray) the C-tail of GRK2, an effective Gβγ scavenger. On the day of the experiment, cells were exposed to the indicated ligands (10 μM; 8 h) or vehicle (DMSO, 0.1% (v/v)). Following washout of the treatment drugs, the cells were exposed to forskolin (6.3 μM) to estimate the sensitization of the cAMP response. Results correspond to the mean ± SEM of the forskolin response in vehicle treated cells transfected with pcDNA3 (B). Cells expressing DOPr and the BRET biosensor as in (A) were treated with pertussis toxin (PTX: 100 ng/mL; 16 h) prior to completing the experiment as in (B). Results correspond to the mean ± SEM of the forskolin response in transfected vehicle treated cells that were not exposed to PTX (C). The number of independent experiments per condition are indicated at the bottom of each histogram bar. Statistical comparisons in (B) and (C) were done by means of one-way ANOVA followed by Sidak’s multiple comparison’s posthoc test, and the results are shown in the figure.

Figure 7

Met-ENK and TIPP induce sensitization of the cAMP pathway via distinct signals. Forskolin concentration response curves show forskolin-induced cAMP accumulation in cells exposed to Met-ENK (A–D), TIPP (E–H) (10 μM; 8 h), or vehicle (CTL) in the presence (white) or absence (gray) of indicated pathway blockers. The results for the concentration response curves represent the mean ± SEM normalized to the maximal forskolin response in cells pre-exposed to vehicle (DMSO, 0.1% (v/v)) in the absence of blocker (CTL). Insets: AUC ± 95% CI values derived from the corresponding curves. The number of independent experiments per condition are indicated at the bottom of the corresponding histogram bars. Statistical comparisons between AUC values obtained in the presence and absence of blocker were computed by verifying overlaps of confidence intervals: CI95 overlaps (no star), CI99 overlaps but not CI95 (*), CI99.9 overlaps but not other CI values (**), and CI99.9 does not overlap (***).

Taken together, the series of results presented above confirm the association between Gα/Gβγ signaling similarities and cyclase superactivation by Ca2+-mobilizing ligands. The validity of this association initially established via signaling similarities is further supported by the observed correlation between AUCs for forskolin-induced cAMP production by Ca2+-mobilizing agonists and AUCs of concentration response curves describing Ca2+ mobilization by these same agonists (Table S4). In contrast, when we tried to correlate superactivation to actual parameters derived from Ca2+ concentration response curves (Emax, Log(τ), or Log(τ/KA)), no significant correlation was detected (Table S4). While numerical imprecisions related to curve fitting could be at the basis of this incongruence, it is important to note that similarities among curve parameters overcame this possible limitation and allowed one to successfully associate specific signals that drive cyclase superactivation to the magnitude of the adaptation induced by Ca2+-mobilizing ligands.

PKC can enhance cellular production of cAMP either via direct phosphorylation of ACs II, IV, V, and VII[64−67] or via the stimulation of the serine threonine kinase Raf-1, which in turn promotes superactivation of ACs V/VI, a mechanism already described for DOPr agonists.[68−70] ACV and ACVI have been found in HEK293 cells[60] (Figure S3). Hence, we determined if the Raf-1 inhibitor GW5074[71] (1 μM) interfered with the sensitization of cAMP production by DOPr agonists. GW5074 reduced Met-ENK-induced sensitization by almost 50% (Figure C) and practically abolished sensitization by TIPP (Figure G), pointing to Raf-1 as a convergence effector for cyclase modulation by Ca2+ mobilizing and nonmobilizing ligands. The fact that superactivation by TIPP relied on Raf-1 but not on the Ca2+/calmodulin/PKC cascade implies alternative mechanisms by which this DOPr agonist may engage Raf-1 to promote cyclase adaptations. In this sense, it is worth considering that TIPP activates the nonreceptor tyrosine kinase Src to produce prolonged ERK activation[19] and that Src also activates Raf-1.[70,72] On the basis of this knowledge, we tested the Src inhibitor PP2 (10 μM) to assess if Src contributed to superactivation by Met-ENK and TIPP. PP2 abolished superactivation by TIPP (Figure H), confirming that Src and Raf-1 underlie adaptations by this agonist. In contrast, in cells treated with Met-ENK, the Src inhibitor PP2 enhanced rather than blocked cAMP sensitization. The contrasting effect of PP2 on the modulation of Src activity by high and low efficacy DOPr agonists has been previously reported and represents the differential ability of high and low efficacy ligands to engage Src-dependent desensitization of the receptor.[19,73] In particular, while TIPP induces long lasting activation of Src and of its downstream pathways, more efficacious DOPr agonists quickly desensitize Src signaling.[19,73] Like TIPP, morphine behaves as a weak DOPr agonist[36] and produces sustained Src-dependent signaling.[19] In keeping with these signaling similarities, the superactivation of cAMP production by morphine was comparable to that produced by TIPP (Figure S5).

Taken together, the results above indicate that Raf-1 is a common effector in the sensitization of cAMP production by the prototypical Ca2+-mobilizing agonist Met-ENK and nonmobilizing agonist TIPP. Through stimulation of Gαi/o, Met-ENK induces the release of Gβγ, the mobilization of Ca2+, and activation of PKC. Via PKC, Met-ENK may drive Raf-1 activity and subsequently modulate ACV and VI.[71−73] Ca2+ mobilized by Met-ENK equally recruits calmodulin, which enhances the activity of Ca2+-sensitive cyclases.[59,64−66] TIPP fails to induce the Ca2+-dependent portion of the sensitization response displayed by Met-ENK, but its weak partial efficacy allows TIPP to elude Src-dependent desensitization of DOPr signaling,[19,73] leading to superactivation of the cyclase path via a mechanism that relies on Gαi/o, Src, and Raf-1. Unlike TIPP but similar to other efficacious DOPr agonists,[19] Met-ENK does not avoid Src-dependent desensitization of the receptor, and this mitigates its overall sensitization of the cAMP pathway. Figure S6 summarizes signals involved in the sensitization of cAMP production by Met-ENK and TIPP.

Ca2+-sensitive cyclases ACI/III[57,59,64−66] are expressed in HEK293 cells[60] (Figure S3) and likely to mediate the observed adaptations driven by the prototypical Ca2+-mobilizing agonist Met-ENK. ACV and -VI that are also expressed in HEK293 cells[60] (Figure S3) are known targets of Raf-1 and are likely crucial for sensitization by the weak partial agonist TIPP. The direct translation of the results obtained in HEK293 cells to in vivo development of tolerance will depend on the expression of the diverse components involved in the supersensitization process in the actual target cells. In this sense, it is worth considering that Ca2+-sensitive cyclases in the central nervous system actively contribute to opioid analgesic tolerance.[29] On the basis of this observation, we expect Gα/Gβγ signaling similarities among Ca2+ mobilizing agonists to have translational value as predictors of cyclase superactivation in target cells expressing ACI, ACIII, and/or ACVIII. Alternatively, when more pertinent cellular models become available for routine drug screening, the same classification process presented herein can be directly applied to the actual target cells of interest.

Materials and Methods

Opioid Ligands

l-Tyrosyl-d-alanyl-l-phenylalanyl-l-alpha-glutamyl-l-valyl-l-valyl-glycinamide (Deltorphin II) was from AnaSpec, and (N-(l-tyrosyl)-(3S)-1,2,3,4-tetrahydroisoquinoline-3-carbonyl)-l-phenylalanyl-l-phenylalanine (TIPP) was from Cedarlane; l-tyrosyl-glycyl-glycyl-l-phenylalanyl-l-methionine (Met-enkephalin), (d-penicillamine2,5)-enkephalin (DPDPE), and N,N-diethyl-4-(phenylpiperidin-4-ylidenemethyl) benzamide (AR-M1000390) were purchased from Sigma-Aldrich. ((+)-4-[(alpha-R)-alpha-((2S,5R)-4-Allyl-2,5-dimethyl-1-piperazinyl)-3-methoxybenzyl]-N,N-diethyl-benzamide) (SNC-80) was obtained from Tocris Cookson.

Chemicals and Reagents

The following chemicals and reagents were purchased from Sigma-Aldrich (St. Louis, MO, USA): pertussis toxin (PTX), 2,2′-(ethylenedioxy) dianiline-N,N,N′,N′-tetraacetic acid (BAPTA-AM), 3-(3,5-dibromo-4-hydroxybenzyliden)-5-iodo-1,3-dihydroindol-2-one (GW5074), 1,2-dimethoxy-12-methyl-[1,3]benzodioxolo[5,6-c]phenanthridin-12-ium (chelerythrine), and 1-[bis(4-chlorophenyl)methyl]-3-[2-(2,4-dichlorophenyl)-2-[(2,4-dichlorophenyl)methoxy]ethyl]imidazol-3-ium (calmidazolium). 4-Amino-5-(4-chlorophenyl)-7-(t-butyl) pyrazolo[3,4-d] pyrimidine (PP2) was acquired from Calbiochem.

Plasmids and Constructs

Receptor Constructs

We used previously designed pSig-Flag-DOPr and pSig-Flag-DOPr-GFP10 constructs.[36] Briefly, pSig-Flag-DOPr was generated by the addition of a signal peptide (pSig) of influenza hemagglutinin (MKTIIALSYIFCLVFA) and the Flag epitope (MDYKDDDDA) at the N-terminal domain sequence of rat DOPr. GFP10 was subcloned in the frame to the C-terminus of pSig-Flag-DOPr to generate the pSig-Flag-DOPr-GFP10 tagged receptor. All constructs were subcloned into pLVX-IRES-Puro.

cAMP Biosensor

The GFP10-Epac-RlucII BRET2-cAMP biosensor[74] is henceforth referred to as the EPAC biosensor. Briefly, the EPAC biosensor consists of human EPAC1 (residues 144–881), which mutated (T781A and F782A). The amino-terminal and carboxy-terminal of human Epac1 are joined by a 5 amino acid residue (GSAGT) linker to Green Fluorescent Protein10 (GFP10) and a 5 amino acid residue linker (KLPAT) to RlucII, respectively. The whole construct is inserted in pcDNA3.1/Zeo (Invitrogen).

Ca2+ Biosensor

Ca2+-regulated photoprotein obelin was cloned into pcDNA3.1/zeo(+).[75]

G Protein Biosensors

Gαi1-RlucII and GαoA-RlucII plasmids were previously described by Gales et al.[37] and by Richard-Lalonde et al.,[39] respectively. Gαi2-RlucII, Gαz-RlucII, Gγ2-RlucII, and Gγ2-GFP10 were described in https://patents.justia.com/patent/9029097. Gβ1 was purchased from Missouri University of Science and Technology (cdna.org). Constructs encoding RlucII-βarr1 and RlucII-βarr2 were, respectively, described by Zimmerman et al.[76] and Quoyer et al.[77] cDNA clones for the following constructs were generously provided as follows: GRK2 and GRK6 by Dr. Antonio De Blasi (Istituto Neurologico Mediterraneo Neuromed, Pozzilli, Italy); GRK5 by Dr. Robert Lefkowitz (Duke University, Durham, NC); Kir3.1 subunit and Kir3.2-GFP10 by Dr. Terry Hebert (McGill University, Montréal, Canada). All constructs were confirmed by DNA sequencing.

Cell Lines and Transfections

HEK293 cells were cultured in 100 mm Petri dishes (Sarstedt, Germany) at 37 °C and 5% CO2 in the Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 10% fetal bovine serum, 2 mM l-glutamine, and 100 unit mL–1 penicillin–streptomycin. For transient transfections of DOPr- and BRET-based biosensors constructs, HEK293 cells were seeded at 3–3.5 × 106 cells/100 mm Petri dish and were grown for 18–24 h before transfecting with polyethylenimine (PEI) (Polysciences Inc., Warrington, PA, USA) at a 3:1 PEI/DNA ratio as per the manufacturer’s instructions. Monoclonal cell lines stably expressing DOPr and the EPAC biosensor (hereafter referred to as EPAC DOPr HEK293 cells) were established by first transfecting 6 μg of the pSig-Flag-DOPr DNA construct/100 mm Petri dish using Lipofectamine (Invitrogen), followed by a puromycin selection (1 μg/mL). This stable cell line was subsequently transfected with 3 μg of the EPAC biosensor, using Lipofectamine (Invitrogen) for the transfection and hygromycin (50 μg/mL) for selection.

BRET Assays

cAMP Accumulation Assays

48 h before the assay, EPAC/DOPr HEK293 cells were plated in 96-well plates at a density of 30 000 cells/well. On the day of the experiment, cells were changed into Tyrode’s buffer (140 mM NaCl, 2.7 mM KCl, 1 mM CaCl2, 12 mM NaHCO3, 5.6 mM d-glucose, 0.5 mM MgCl2, 0.37 mM NaH2PO4, 25 mM HEPES, pH 7.4) and incubated for 30 min at 37 °C in 5% CO2 before starting the manipulations. For experiments designed to evaluate how prolonged exposure to DOPr agonists influenced cellular production of cAMP (superactivation assays), EPAC/DOPr HEK293 cells were incubated for 8 h in the presence of Met-ENK, TIPP, deltorphin II, DPDPE, ARM390, SNC-80 (10 μM), or vehicle (DMSO; 0.1% (v/v)). This treatment concentration ensured that the superactivation responses were obtained at a maximal effective concentration by all ligands. The treatment duration was established in pilot studies where 4 h of treatment did not have an effect while 8 h of exposure resulted in clear, measurable differences in sensitization induced by different agonists. At the end of the treatment, the cells were washed with Tyrode’s buffer (3× for 3 min at 37 °C) and were then redistributed into 96-well plates (PerkinElmer, Waltham, MA; 2 × 104 cells/well). Cells were then incubated for 8 min with coelenterazine 400a (5 μM) (Bioshop, Canada) and forskolin (Bioshop, Canada) before taking BRET measures at 37 °C. We have previously shown that within this delay the BRET signal monitoring cAMP levels reaches a plateau.[38] BRET2 signals were determined by calculating the ratio of the emission at 530 nm (GFP10) over the emission at 400 nm (RlucII) using a Mithras LB 940 Multimode Microplate. When pathway inhibitors were used to block the superactivation of cAMP by the DOPr agonists, these were introduced throughout exposure to the agonist. Preliminary experiments were carried out to determine the minimal concentrations needed to produce the maximal inhibition of the cAMP response (data not shown). The results showed that the concentration of inhibitors used were as follows: BAPTA-AM (3 μM), calmidazolium (10 μM), chelerythrine, (5 μM), PP2 (10 μM), and GW5074 (1 μM). When assessing efficacy of different ligands to modulate cAMP accumulation, coelenterazine 400a was added to the cells as above (5 μM, 3 min), followed by forskolin (15 μM, 3.5 min) and increasing concentrations of DOPr agonists. BRET2 readings were taken as above, 5 min after ligands were introduced.

G Protein Activation Assays

HEK293 cells cotransfected with rat DOR, Gβ1, and Gγ2-GFP10 and of each of the different Gα subunits (Gαi1, Gαi2, GαoA, and Gαz) tagged with RlucII were incubated with coelenterazine 400a (5 μM for 3 min) before exposing them to increasing concentrations of DOPr agonists for an additional 5 min. BRET readings were taken as above.

β-Arrestin2 Recruitment Assays

Cells were cotransfected with rat DOPr-GFP10 and either β-arrestin1-RlucII or β-arrestin2-RlucII with or without human GRK2, GRK5, or GRK6. Cells were stimulated for 10 min with increasing concentrations of DOR ligands and then incubated with coelenterazine 400a (2.5 μM for 5 min). β-Arrestin2 recruitment was assessed using BRET2 filters.

Kir 3.2 Channel Activation Assay

HEK293 cells were cotransfected with constructs encoding rat DOR, human Kir3.2-GFP10, and human Gγ2-LucII. Cells were stimulated for 5 min with increasing concentrations of DOR ligands and then incubated with coelenterazine 400a (2.5 μM for 5 min). Channel activation was measured in BRET2.

Ca2+ Mobilization

As described previously,[78] HEK293 cells were cotransfected with rat DOR and obelin. On the day of the experiment, the cells were washed and preincubated with coelenterazine cp (1 μM) and kept in the dark at 25 °C for 2 h. DOR ligands were subsequently injected into the wells, and the blue luminescence emission (465–495 nm) was monitored every 0.5 s for 60 s using a SpectraMax L microplate reader (Molecular Devices, Sunnyvale, CA). Results were expressed as relative luminescence units (RLUs). The AUC of 60 s stimulation by agonists was calculated and normalized to Met-ENK’s maximal response.

Data Analysis

Superactivation of the Adenylyl Cyclase Pathway

The superactivation of the adenylyl cyclase pathway was evaluated as the change in forskolin-induced cAMP production following the exposure of cells to the vehicle, Met-ENK, or other the indicated DOPr agonists. Typical independent experiments comprising concentration response curves for forskolin were performed in duplicate for each condition. For each independent experiment (one of N), we calculated the average of the replicates and normalized responses across the experiments, fitting to the logistic equation as follows. All data were fit simultaneously, and concentration response curves from each independent experiment were fit to the equation, where EC50 was fit as a logarithmand the following constraints were applied: (i) each curve from the same drug had the same span (max asymptote minus min asymptote), EC50, and Hill coefficient and (ii) each curve from the same experiment shared the same normalization factor (panel scale, which was applied to the span). During the minimization (fitting), we computed AUC as a transformed parameter (a parameter that is computed from the other parameters and not from the data itself). At each iteration, the AUC for the estimates was computed using numerical integration (c language routines from Press et al.[79]) in the concentration range of 10–9–10–4 M (but on the common Log10 scale in order for AUC values to be meaningful). In effect, this approach outputs an estimate of the variance of the AUC, which as for the other parameters, can be used to compute its standard error.

Estimates for confidence intervals at 95% confidence levels (CI95) (also at CI99 and CI99.9) and standard errors (SE) were computed as described by Motulsky and Christopoulos.[80] For each estimate, the variance is used to compute the SE from the squared sum of residuals (SSR) and the degree of freedom (df).For each parameter, there are corresponding specific values of variance, SSR, and df. An SSR corresponds to each concentration response curve. The SSR for parameters estimated from multiple curves is their sums, and their dfs are the count of data points in them minus the number of estimates minimized using them. t tests can be conducted for each estimate using these SE and df values. The values of the confidence intervals also are computed from the SE and df as follows. Confidence interval values are obtained by multiplying the SE by the Student’s t value corresponding to df and alpha (where alpha = 0.025 for CI95, 0.005 for CI99, and 0.0005 for CI99.9).

Curve Fitting for Signaling Profile

Concentration response curves for signals monitored with different BRET biosensors were fit with both the logistic equation (as described above) and the operational model.[81]

The way the operational model can be applied to assess functional selectivity of the GPCR ligands is being actively pursued, particularly concerning the use of affinity information for minimizing concentration response curves obtained in cell-based bioassays.[82−85] One of the proposed uses of the model posits that the receptor fully uncoupled from downstream transducers/effectors supports all signals generated by a given ligand–receptor pair. On the basis of this reasoning, the fitting method estimates Log(τ) values by constraining KA to experimentally obtained affinities.[86] A limitation to this use of the model is the confidence with which experimentally determined affinity values describe the interaction between the ligand and the fully uncoupled form of the receptor.[84] An alternative and more frequent use of the model adopts no constraint on KA values, except in the case of full agonists where very high affinity is assumed.[82,87,88] A consequence of this minimizing strategy is that Log(τ) (and KA) values for full agonists remain undetermined such that these ligands are solely described in terms of transduction coefficients (Log(τ/KA)).[82,83] Not being able to obtain efficacy (Log(τ)) estimates for full agonists poses a problem in view of classifying ligands according to signaling similarities across comprehensive signaling profiles.[36] To circumvent this problem, we reasoned that curve fitting of fully effective responses could be guided by considering KA information from all of the ligand’s signaling readouts, including those where it displays partial efficacy. This reasoning was embodied by constraining KA values to be shared across all functional readouts for each ligand–receptor pair. For this purpose, operational parameters for all curves of a specific drug–receptor pair were minimized in a single execution so that a best estimate of KA, which simultaneously satisfied best fits for all curves, was obtained. At the same time, the other parameters (including Log(τ)) were minimized for curves describing each specific functional readout, unhindered by sharing of KA. The fitting procedure is similar to the one described for shared KA values across curves obtained at varying receptor densities,[81,89] and the end result is an estimation of sensor-specific efficacy Log(τ) and ligand–receptor affinity (KA) values.

Importantly, the affinity estimates that are obtained with a shared KA across functional responses are not equivalent to the “functional affinities” that are estimated by the more frequent unconstrained use of the model.[82,83] However, both types of estimates were reasonably correlated, as verified by comparing shared KA values to those obtained with the unconstrained approach (r2 = 0.7736; p = 0.02).[83] Furthermore, since neither “functional affinities”[83] nor shared KA estimates necessarily represent actual binding parameters to the uncoupled form of the receptor, Log(τ) values obtained with either of these methods are not equivalent to the operational efficacy estimates obtained using experimental binding data to assign KA values.[82,83,86] A limitation of the approach we propose concerns the minimization of the curves for highly efficacious agonists that do not display a partial response at any of the readouts tested. In such cases, the indetermination for KA and Log(τ) values would persist across all readouts. For the group of ligands in the current study, Log(τ) (and KA) estimates for all ligands were obtained though the number of partial responses available for each ligand determined the fitting error of these parameters (see Table S1). In summary, shared KA minimization allowed us to obtain information to classify ligands according to Log(τ) information without introducing the affinity confounder present in Log(τ/KA) coefficients.

Clustering Drugs According to Signaling Profiles

Drugs were grouped according to the similarities in parameters generated from their concentration response curves (Figure A–C). To do so, we used a previously described statistical method whose output is a similarity matrix describing the frequency of coclustering of pairs of ligands across iterative comparisons of parameters built into the procedure.[36] Briefly, the method performs as follows. Given the matrices of fit parameters and the corresponding error estimates, it generates a number of replicate matrices by sampling the underlying distribution of fit parameters and subsequently submits each of these matrices to a NMF factorization followed by multiple K-means clustering. The K-means clustering frequencies are interpreted as similarity measures, which are averaged. Finally, using R, the average similarity matrices containing the frequency of the coclustering values are used as input to the heatmap function. This last step computes hierarchical clustering of the drugs whose output tree is shown alongside the reordered similarity matrix as an intuitive heatmap. The resulting pairs of trees/heatmaps are shown in Figure . By taking the row corresponding to the reference compound (Met-ENK), we obtain the frequency of coclustering of every compound in relation to it. This reveals how similar each compound is relative to the reference compound. The jupyter notebook python script IterativeClustering_NMF.ipynb, which is part of the github package https://github.com/JonathanGallion/Benredjem-Gallion, was used to compute similarity matrices using 100 replicates and 25 NMF restarts (numits).