BiasNet: A Model to Predict Ligand Bias Toward GPCR Signaling

Model Performance.

A total of five different machine learning algorithms along with graph-based neural networks (see Machine Learning section) were trained on the BiasDB data using 15 different descriptors. The best-performing models based on F1 score are reported. The test set scores for each of our five algorithms are shown in Figure 1. For the six metrics used to evaluate our models (F1 score, Recall, Precision, Specificity, Accuracy, and AUC), scores that are closer to 1 indicate better models. Although there is no strict threshold for a “good” model, we can reasonably assume our models were able to correctly classify BiasDB ligands as G protein biased or β-arrestin biased because all metrics for the top-performing models (with the exception of D-MPNN) were above 0.70. All of our models scored an accuracy and AUC above 0.8. Also, all of our models showed a recall above 0.70 and a precision above 0.65. The F1 score for all of our models was above 0.75. The specificity of our models was above 0.8, and for two of our models (XGB and RF), this number was above 0.9. McNemar test results (Table 1) show there is no statistically significant difference between our five models (all pairwise p-values are >0.05). Based on these results, we conclude our models were able to distinguish β-arrestin biased ligands from G protein biased ligands in our training set. Additionally, we observe there was no statistically significant difference between our models’ abilities to distinguish G protein biased ligands from β-arrestin biased ligands.

With these results in mind, it is important to recognize that G protein/β-arrestin bias is a continuum. While the molecules in BiasDB are grouped together in convenient categories that lend themselves to train a classification model (G protein or β-arrestin), it is possible that similar regression models which predict the level of bias can be built. We did not pursue this goal primarily because the metrics used to quantify bias in BiasDB are not consistent from ligand to ligand. For each ligand, two assays are required to determine bias: one that establishes β-arrestin activity and one that establishes G protein activity. While the β-arrestin recruitment assay can be used consistently to determine β-arrestin bias, many different assays can be used to gauge G protein activity (cAMP production, GTPase activity, BRET, and various others). It is because of this reason that we restricted the scope of this study to the classification problem.

External Predictivity.

Having discussed the reliability of our models on BiasDB, we also wanted to demonstrate comparable performance with external test data, thus demonstrating the generalizability of our models. We considered compounds that show potentially therapeutic GPCR signaling bias.³ Out of the compounds discussed in the work done by Tan and colleagues, 21 were not found on BiasDB. The bias of these 21 compounds was predicted using our models. The list of compounds is shown in the Supporting Information. The bias of 18 of these compounds—approximately 90%—was correctly predicted. This finding corroborates the performance of our models on data not found in BiasDB.

Applicability Domain.

Determining whether a compound is biased toward G proteins or β-arrestin is only appropriate when the compound targets GPCRs. To address this concern, the web implementation of our models provides the user information about whether the compound is likely to bind to a GPCR. We developed an SVM machine learning model trained on GPCR-interacting and non-GPCR-interacting ligands to restrict the applicability domain to cases where G protein bias is relevant. Model results are shown in Table 2. Here, we note the applicability domain model used by BiasDB is intended primarily to serve as a warning for the user. The GPCR-binding model is trained on ligands, which bind at the orthosteric sites as well as allosteric sites.

Molecular Fragments Related To Bias.

One major goal of our work was to understand intuitively what might make a molecule biased. Molecular fingerprints encode the chemical groups of a molecule into a vector and aid in training accurate models; however, the meaning of this vector is not readily understood by humans. To address this issue, we developed an MLP model trained on RDKit descriptors. The top-20 most important RDKit descriptors were analyzed, and from this list, we observed the two most important fragments are the number of secondary amines and the number of aromatic amines.

Descriptive statistics summarizing the number of secondary amines and aromatic amines for BiasNet ligands are shown in Table 3. The table shows ligands in BiasDB had 0 to 14 secondary amines. While this first range is large, it is worth noting that 75% of G protein biased ligands had zero or one secondary amine and 75% of β-arrestin biased ligands had zero, one, or two secondary amines. The range of aromatic amines is markedly smaller, however. This range is 0–3 aromatic amines. We observe aromatic amines are also more common in β-arrestin biased ligands. In particular, the value at Quartile 3 (Q3) for G protein biased aromatic amines is 0, and the value at Q3 for β-arrestin aromatic amines is 1.

The mean number of secondary amines for G protein biased ligands was 1.07, while for β-arrestin biased ligands, it was 2.44. The mean for G protein aromatic amines was 0.22, and for β-arrestin, this number was 0.49. One-tailed Z-tests comparing the difference in means for secondary amines and for aromatic amines show a statistically significant difference between the averages of β-arrestin biased ligands and G protein biased ligands (p-value < 0.01 for both of these fragments).

Figure 2 more completely shows the distribution of molecular fragments. The violin plot on the left of Figure 2a exhibits one high-density region at the base, corresponding to 0–1 secondary amines. The violin plot on the right has two such regions. The first region is at the base of the plot and has the greatest width at around 1 secondary amine. The second region is at the top of the plot and has the greatest width at 10 secondary amines. Figure 2b shows that G protein biased ligands had the greatest density at 0 for aromatic amines. Similarly, β-arrestin biased ligands had the greatest density at 0 aromatic amines. The difference is that for β-arrestin biased ligands, a greater proportion of compounds have 1 or 2 aromatic amines.

Taken together, these statistics suggest increased secondary amines and aromatic amines might be connected to β-arrestin bias relative to G protein bias. While our models did identify this pattern, we recognize that other features play a larger role in determining bias. Thus, it is not possible to say with complete certainty that a molecule with no secondary amines or no aromatic amines is G protein biased. Similarly, it is not possible to say that a molecule with many secondary amines or aromatic amines is β-arrestin biased—merely that these fragments tend to be found in greater abundance in β-arrestin biased ligands.

In addition to highlighting these statistics, our aim was to try to understand the structural reasons that might lead to bias. We focus our discussion on secondary amines and consider the two molecules in Figure 3. Both molecules are ethanolamines and agonists for the β-2 adrenergic receptor (β2R), and have a similar mechanism of action. The molecule on the top panels (isoproterenol) contains only one secondary amine, while the molecule on the bottom panels (BI-167107) contains two secondary amines. The molecule on the top panels is unbiased, while the molecule on the bottom panels is β-arrestin biased.

Previous docking studies have demonstrated that isoproterenol can form hydrogen bonds between its para hydroxyl and SER 203, SER 204, and SER 207 on Transmembrane 5 (TM5) of the β-2 adrenergic receptor.⁸ Additionally, the meta-hydroxyl group can form hydrogen bonds with SER 203. We suspect weak interactions with TM5 might be connected to β-arrestin bias.

We performed docking studies on isoproterenol and compared them to BI-16707 co-crystallized to the β-2 adrenergic receptor. Out of all docking poses, we chose the one with the second-lowest binding energy as it closely resembles the pose reported in the literature by Weiss and colleagues.⁸ All binding poses and associated binding energies are available in the Supporting Information Figure S1. The binding poses of the two compounds are shown in Figure 3. Like in the study by Weiss et al., the para-hydroxyl group on isoproterenol appears in close proximity to SER 203 and SER 204; however, our model does not find a hydrogen bond between the para-hydroxyl group and SER 204. Still, our simulation did find a hydrogen bond between the meta-hydroxyl and SER 203. Importantly, in both the docking simulations done by Weiss et al.⁸ and in this study, isoproterenol explicitly interacts with TM5.

Co-crystallized BI-16707 shows similar proximity to SER 203 and SER 204 with its para-hydroxyl group and a secondary ring member amine; however, the co-crystallized structure did not demonstrate any hydrogen bond. While it is still possible that BI-16707 interacts with TM5 residues via dipole–dipole interactions, these are not as strong as hydrogen bonds.

Another factor to consider is the role the ring member nitrogen plays in determining bias. Whether this nitrogen participates in dipole–dipole interactions or hydrogen bonds with residues from TM5, such interactions would be weaker than analogous interactions with a hydroxyl from the isoproterenol group due to the higher electronegativity of oxygen.

We believe the stronger interaction of isoproterenol with TM5 might account for the balanced bias demonstrated by isoproterenol. By the same token, the weaker interactions between BI-167107 might explain β-arrestin bias. Experiments showing isoproterenol co-crystallized with the β-2 adrenergic receptor are needed to clarify the role of SER 203, SER 204, and SER 207 and confirm our hypothesis about β-arrestin bias.

Unsupervised Learning.

To corroborate the findings of our supervised machine learning models, we performed an unsupervised machine learning technique, t-stochastic neighbor embedding (t-SNE). We consider the 5 clusters in our discussion represented in Figure 5. Templates for each of the five clusters are shown in Figure 6. A full collection of cluster compounds can be found in the Supporting Information.

Cluster 1 (enclosed by a red box in Figure 5) identifies 15 G protein selective EP2 agonists. These ligands (represented in Figure 6) share three features that form a common structural backbone. The SAR of these compounds has previously been described by researchers synthesizing novel EP2 selective agonists.¹¹ The first feature is a thiazole ring, which is important for EP2 receptor subtype selectivity.¹² The second feature is a middle bicyclic moiety, which allows the molecules to mimic endogenous prostacyclins. The third feature is a variable ω chain (corresponding to the R1 group), which contributes to bias. Importantly, we note that these ligands do not contain secondary amines consistent with the idea that G protein biased ligands tend to have less of these moieties than β-arrestin biased compounds.

Cluster 2 identifies G protein biased κ opioid receptor agonists. These ligands share a common central triazole which forms a scaffold attached to three different aromatic rings. One of these rings is a 4-pyridyl ring, and the remaining two are variable aromatic rings. The SAR of these compounds has previously been described by researchers investigating biased κ opioid receptor agonists.^13,14 To alter bias, the non-pyridyl rings are exchanged, but the core triazole ring is left intact. The core triazole also contains no secondary amines.

Cluster 3 represents G protein biased D1 receptor agonists. Of the 12 ligands in this cluster, all are G protein biased. Researchers have previously described the binding mode of these compounds, referred to in the literature as non-catechol D1 agonists.¹⁵ These ligands interact with SER 188 and LEU 190 from intracellular loop 2 (ECL2) via the variable R3 group.

Cluster 4 represents β-arrestin biased D2 receptor agonists. Cluster 4 contains 17 ligands, 16 of which are β-arrestin biased. Aripiprazole is the exception to this cluster and is the unbiased reference compound by which biased compounds are compared. The SAR of these compounds has previously been characterized.^10,16–18 Each of these ligands can be divided into four regions: a “left-hand side (LHS) phenyl ring,” “middle cyclic amino moiety,” “central linker,” and “right-hand side (RHS) bicyclic aromatic moiety.”¹⁰ The region we will discuss is the RHS bicyclic aromatic moiety pictured in Figure 6. This moiety contains one secondary amine. Interestingly, the work done by Chen and colleagues found one compound (not included in this cluster) for which the RHS bicyclic aromatic moiety is replaced with a benzimidazolone group. This compound, 5-(4-(4-(2,3-dichlorophenyl)piperazin-1-yl)-butoxy)-1H-benzo[d]imidazol-2(3H)-one, has two secondary amines and shows an increased β-arrestin bias relative to aripiprazole, further highlighting the connection between secondary amines and β-arrestin bias.

Cluster 5 contains 19 ligands, 18 of which are G protein biased. Among the receptors targeted by cluster 5 are the μ-opioid, κ-opioid, nociceptin, apelin, and ghrelin receptors. We limit our discussion to μ-opioid receptors agonists, as the majority (13) of ligands from cluster 5 are μ-opioid agonists. The μ-opioid ligands in cluster 5 are G protein biased and have the common structure, as shown in Figure 6. All variations of these compounds have at most 1 secondary amine—the amine in the 5-member ring. While this finding does not definitively require that all G protein biased μ-opioid receptor targeting compounds have at most 1 secondary amine, it further highlights the connection we observed between a small number of secondary amines and G protein bias.

Clinically Significant Compounds with HierS Scaffolds.

To end our discussion of scaffolds, we chose to illustrate how bias can have a clinically significant effect. Here, we consider three commercially available medications, as shown in Figure 7. Oliceridine is a recently FDA-approved G protein biased drug molecule containing scaffold 1. Oliceridine is a μ-opioid receptor agonist intended to treat severe acute pain. β-Arrestin activity for μ-opioid receptors is associated with increased gastrointestinal dysfunction, respiratory dysfunction, and decreased efficacy.¹⁹ Because oliceridine shows a marked G protein bias, researchers suspect oliceridine might have a better safety profile than commercially available drugs like morphine. Despite this optimism, other researchers suggest the enhanced safety profile of oliceridine might be attributed to a lower intrinsic activity for the drug.²⁰ In any case, more research is needed to confirm whether oliceridine is a suitable alternative for current opioid drugs.

Lisuride is a β-arrestin biased drug containing scaffold 3. Lisuride binds to D2 receptors and is indicated for Parkinson’s disease. Standard treatments for Parkinson’s include the drug l-Dopa, which activates dopamine receptors in the body. l-Dopa is associated with side effects such as dyskenesia, which are thought to be caused by receptor supersensitivity.²¹ Lisuride bypasses such side effects because its β-arrestin action mitigates this receptor supersensitivity. For this reason, lisuride is often prescribed before l-Dopa.

Clonidine is a β-arrestin biased drug containing scaffold 4. Clonidine binds to the α2 adrenergic receptor and is indicated for high blood pressure. Work has been conducted for related groups of anti-hypertensive and cardio-protective drugs targeting the β1 adrenergic receptors. This work links β-arrestin activation to “therapeutically advantageous” activity.²² Although clonidine works through a different receptor and overall mechanism of action, there is reason to believe β-arrestin bias for clonidine could also be a desirable characteristic for the drug. For clonidine, β-arrestin bias is liked to increase endocytosis of the α2 adrenergic receptor.²³ This β-arrestin-directed effect presumably leads to desensitization of the ligand, which results in decreased drug side effects.

Molecular Features.

Supervised machine learning models included those trained solely on molecular descriptors and one trained solely on scaffolds. In total, 15 molecular descriptors (each represented by a vector of different size) were used to train each of 15 descriptor models, and one vector with 704 scaffolds was used to train a scaffold model. The model used to make predictions on the Biasnet web service uses the long ECFP4 vector described below. The scaffolds used to train the scaffolds model are discussed in the following Scaffolds section. The unsupervised t-SNE model was trained solely on RDK7 fingerprints, as the t-SNE algorithm relies on dimensionality reduction, and using multiple descriptors would make this task more difficult. Here, we discuss each of the descriptors in more detail.

Scaffolds.

All queried G protein biased ligands and β-arrestin biased ligands were introduced into the hscaf web app available at http://pasilla.health.unm.edu/tomcat/hscaf/hscaf. After processing, 704 scaffolds were generated. Thus, a 704-bit vector was used to train a RF model, where the presence of a scaffold was denoted by a 1 in the vector, and the absence was denoted by a 0. The scaffolds are included in the Supporting Information.

Following model development, importance scores were calculated using the permutation feature importance module found in Scikit-learn. Scaffolds with a score greater than 0.005 were chosen for screening. Selection based on importance score resulted in 10 scaffolds. Out of these, scaffolds were cross-checked against BiasDB to ensure that the scaffolds represented primarily one type of bias. Five scaffolds fit this criterion.

Machine Learning.

The main data set was divided into training (80%) and test (20%) sets using a stratified sampling method, which ensures the equal distribution of different classes into both sets. Firstly, four different machine learning (ML) algorithms, namely random forest (RF), extreme gradient boosting, multilayer perceptron (MLP), and support vector machine (SVM) were trained with their default settings, with 15 different descriptors (see Molecular Descriptors). Finally, the best-suited combination of classifiers and features (RF/Avalon, XGB/Avalon, MLP/lecfp4, and SVM/ecfp2) were selected and further tuned based on hyperparameter optimization. We used Hypeopt,³² where each model was tuned based on the grid search method, with default 3-fold cross-validation. To compare the model performances between the models trained with manually computed molecular fingerprints and models trained on graph-based neural networks that learn by operating on the graph structure of the molecule, we retrained the models using Chemprop D-MPNN. The D-MPNN model was also trained with 3-fold cross-validation. It was optimized using various model hyperparameters such as the hidden size of the neural network layers, models depth, number of feedforward network layers, and dropout probability. Among all of the optimized models, the outperforming model based on the test set F1 score was used for web implementation. The five aforementioned ML algorithms are summarized as follows.

Machine Learning Model Comparison.

Our Machine learning models were evaluated by receiver operation characteristic area under the curve (ROC AUC), accuracy, F1 score, recall, precision, and specificity. Additionally, a pairwise McNemar test was conducted to compare the performance of each machine model against all of the other models. These metrics are calculated with true positive (TP), true negative (TN), false positive (FP), and false negative (FN) classifications.

Applicability Domain.

The target central resource database version 540 (TCRD v540) was used for training the applicability domain model. The compound activity table was merged with the protein table to obtain a mapping of ligand to target. For this combined table, rows with target descriptions containing the phrase “G protein-coupled” were assigned a label of 1. Rows not containing the previous description were given the label 0. The data set was featurized using long ECFP4s. A linear SVM model was trained and tested on these features following a 75/25 train/test split. The regularization hyperparameter C was set to 1 during model development. In total, 381 652 ligands were used, 709 of which showed activity for GPCRs and 380 943 for non-GPCRs.