Abstract
Mathematical models of biochemical reaction networks are an important and emerging tool for the study of cell signaling networks involved in disease processes. One promising potential application of such mathematical models is the study of how disease-causing mutations promote the signaling phenotype that contributes to the disease. It is commonly assumed that one must have a thorough characterization of the network readily available for mathematical modeling to be useful, but we hypothesized that mathematical modeling could be useful when there is incomplete knowledge and that it could be a tool for discovery that opens new areas for further exploration. In the present study, we first develop a mechanistic mathematical model of a G-protein coupled receptor signaling network that is mutated in almost all cases of uveal melanoma and use model-driven explorations to uncover and explore multiple new areas for investigating this disease. Modeling the two major, mutually-exclusive, oncogenic mutations (Gαq/11 and CysLT2R) revealed the potential for previously unknown qualitative differences between seemingly interchangeable disease-promoting mutations, and our experiments confirmed oncogenic CysLT2R was impaired at activating the FAK/YAP/TAZ pathway relative to Gαq/11. This led us to hypothesize that CYSLTR2 mutations in UM must co-occur with other mutations to activate FAK/YAP/TAZ signaling, and our bioinformatic analysis uncovers a role for co-occurring mutations involving the plexin/semaphorin pathway, which has been shown capable of activating this pathway. Overall, this work highlights the power of mechanism-based computational systems biology as a discovery tool that can leverage available information to open new research areas.
Similar content being viewed by others
Introduction
Data intensive computational methods have contributed to significant advances in a variety of scientific disciplines. Cancer biology in particular has seen many advances that follow from the application of computational methods to emerging forms of acquired data1,2,3,4,5,6. However, not all cancers have abundant data. Rare cancers, and rare subtypes of cancer, do not generally have large data sets of the scale required for many computational methods. The number of patient samples, tumor derived cell lines, and mouse models that can be used to generate new, large data sets are commonly also limited for rare cancers. However, years of broad research efforts for individual rare cancers have resulted in several areas where modest amounts of data are available. For example, cell biological, biochemical, genomic (DNA), and transcriptomic (RNA) characterizations are frequently available for rare cancers. We hypothesized that biochemical-mechanism based models of the essential molecular network(s) in a rare cancer may provide a useful approach for re-analyzing the available data to generate new inferences and insights.
Previously, we have used mathematical modeling to investigate signaling by pathogenic RAS pathway mutations7,8,9,10,11,12,13. In these previous studies, we simulated the core biochemical processes that regulate RAS signaling. The model parameterized for wild-type proteins describes the dynamic equilibrium of signaling in a resting cell. Mutant versions of RAS proteins are characterized by changes in these biochemical rate constants and elevated levels of steady-state signaling14,15,16. By substituting biochemical parameters that involve the reactions of a given protein with the biochemical kinetic parameters (e.g., kf, kr, kcat, Km, etc.) the model predicts how signaling outcomes differ when the mutant is present17,18,19. This approach has been useful for uncovering new aspects of RAS mutant biology7,8, for solving long-standing questions about RAS mutant biology10,11,20, and for generating new ideas about targeting the pathway7,9,12,13. This general approach has been reproducibly used by different groups to study oncogenic RAS signaling7,14,16,21,22,23, attesting to the reproducible nature of the general computational approach even if the specific details of the model differ from group to group (just as different experimental groups may utilize slightly different reagents or methods for otherwise conceptually similar experiments). Moreover, the approach of comparing physiological signaling (from wild-type proteins) with pathological signaling (by mutant proteins) has been broadly used for diseases other than cancer and/or for proteins other than RAS24,25,26,27,28,29,30,31,32, further highlighting the potential power of the method to be generalized and applied more extensively. Previous work has also highlighted that mathematical models of a modest scope (approximately five to twenty-five molecular species) can be generated and used to formulate novel hypotheses about pathways important to disease on a time-scale comparable to experimental screens (~6 months)33. Altogether, the previous studies suggest that mathematical models of disease-causing mutations in signaling networks could potentially be useful as a discovery tool that can open new areas of research for subsequent validation of model-based inferences and predictions.
To do this, we combine mechanism-based computational modeling with bioinformatic analysis and experimental biology to investigate uveal melanoma (UM), a form of cancer that has an incidence of approximately five per million individuals in the United States34. We focused on UM as a test case because there are much less data for UM than for more common cancers, because of the importance of G-protein coupled receptor signaling in UM, and because oncogenic G-proteins have proven well-suited for mechanistic modeling. We first developed a biochemical mechanism-based mathematical model of the biomolecular signaling pathway that is mutated and plays a causal role in nearly all cases of UM. As part of our model validation, we demonstrate our model can explain confusing aspects of recently developed direct pharmacological inhibitors against the major driver of this cancer35,36,37,38. Next, we modeled and compared two major classes of UM. Our mechanistic model revealed gaps in available knowledge that are of fundamental importance to understanding the molecular basis of this disease. We addressed these existing gaps in knowledge with new experiments and we updated our model to the new information. These insights led us to hypothesize that there may be another pathway that augments the known pathways driving UM. Evaluation of the available, limited, DNA sequencing data with our focused hypothesis in mind led us to identify statistically enriched mutations in the plexin/semaphorin signaling pathway. Overall, our study demonstrates that this integrated approach, which we call MAGPIE for Model-Assisted Generation of Predictions and Integration of Experimental Data, can catalyze scientific progress through its ability to leverage multiple “modest data” resources.
Results
Development of a mechanistic mathematical model of oncogenic GPCR-pathway signaling in UM
Uveal melanoma (UM) is a rare but deadly cancer of the eye34. Although both UM and cutaneous melanoma arise from melanocytes, the two malignancies are biologically and genetically distinct39. Despite the great progress made in recent years in the development of both targeted therapies and immunotherapies for cutaneous melanoma40, there are no FDA approved targeted therapies for UM39 and only one FDA approved immunotherapy41. The identification of new pathways and targets for the treatment of UM thus remains a priority.
UM is highly dependent upon G-protein coupled receptor (GPCR) signaling (Fig. 1a). GPCRs are one of the largest and most diverse families of membrane receptors in the genome42. GPCRs bind a wide variety of extracellular ligands which result in conformational changes of the receptor. Activated GPCRs catalyze the exchange of GDP for GTP in the alpha subunits (Gα) of heterotrimeric G-proteins downstream. Once bound to GTP, the Gα subunits effectively disassociate from the beta-gamma subunits (Gβγ) to bind and activate effectors downstream. Evaluation of the somatically acquired mutations in UM have highlighted a critical role for one GPCR pathway. Mutually exclusive mutations in the CYSLTR2/GNAQ/GNA11/PLCB4 signaling pathway occur in >95% of UM cases39. Mutations to the Gα encoding genes GNAQ and GNA11 (which will be referred to together as GNAQ/11) account for the overwhelming majority (~90%) of these mutations in UM43.
a Schematic of canonical G protein-coupled receptor (GPCR) signaling. b Schematic of the mechanistic mathematical model of Gαq/11 and CysLT2R signaling in uveal melanoma (UM) with chemical kinetic parameter constants annotated. Within the G protein activation/inactivation cycle, CysLT2R is the guanine nucleotide exchange factor (GEF) that activates the small G-proteins Gαq and Gα11 by promoting the replacement of bound guanosine diphosphate (GDP) for guanosine triphosphate (GTP), regulator of G protein signaling (RGS) is the GTPase-activating protein (GAP) that inactivates the small G-proteins Gαq and Gα11 by promoting the hydrolysis of bound GTP to GDP, and both TRIO and PLCβ are effectors that bind GTP-bound Gαq and Gα11. Within (a, b), the schematized G-protein alpha subunits (red circles) can represent either (or both) wild-type and mutant forms.
The most frequent GNAQ/11 mutations occur in a hotspot at codon 209. Two of these mutations, GNAQ/11 Q209L and Q209P, have been shown to encode mutant Gα subunits that differ in key biochemical properties, including binding to effector proteins and to regulator of G protein signaling (RGS) proteins. RGS proteins are the GTPase-activating proteins (GAPs) that catalyze the conversion of GTP-bound Gαq/11 to GDP-bound Gαq/1144. Although there are significant differences in several of the biochemical properties of WT and Q209L/P mutant Gαq subunits, the importance of these various differences for tumorigenesis in vivo remains unclear. Mathematical modeling of GPCR signaling pathways has a substantial history of aiding in understanding the dynamics of GPCRs in response to ligands and pharmacological targeting45,46,47. We here developed a mechanistic mathematical model of the Gαq/11 pathway that begins at the level of the GPCR CysLT2R (the protein encoded by the CYSLTR2 gene), that includes the G-protein activation cycle, and that includes the binding to its two primary effector proteins, PLCβ and TRIO48,49 (Fig. 1b). The activation of CysLT2R by ligand and the transmission of signals downstream from PLCβ and TRIO are not included explicitly in the model (see below for more details).
Some key features of the proposed model are as follows: first, the reactions are modeled mechanistically using ordinary differential equations (ODEs). This mathematical formalism is based on mass action and classical enzyme kinetics. Direct measurements of kinetic rate constants for the various reactions were taken from the literature wherever possible. Second, a key component of the model is that when considering a particular oncogenic mutation (i.e., the Q209L mutation in Gαq) this mutant protein is considered explicitly as a distinct species in the model with its own biochemical properties. For example, to model the case of a heterozygous GNAQ Q209L mutation, of the total pool of Gα subunits in the model, 75% will be WT and 25% will be mutant (this reflects an assumption of the model that the Gαq and Gα11 subunits are biochemically identical and expressed in similar amounts in vivo, thus Gαq and Gα11 each account for 50% of the total Gα pool in the model). The main biochemical differences between WT and mutant Gαq subunits reported in the literature are the rates of basal and GAP-stimulated GTP hydrolysis and the affinities for effectors and GAPs (the supplement includes a full discussion of the literature concerning biochemical characterization of various mutants). For model parameterization, we initiate our study by identifying a set of approximate values for each parameter based upon the relevant literature (see the supplement for details). Additional mathematical details, references for the biochemical parameters of the model, and a Python notebook that includes all of the code needed to run and reproduce the model are included as supplementary information and/or on github.
While oncogenic activation ultimately depends on downstream targets, such as Yes-associated protein (YAP) for TRIO and the Extracellular Signal-Regulated Kinase (ERK) Mitogen Activated Protein Kinase (MAPK) pathway via protein kinase C (PKC) for PLCβ, the model remains agnostic to these outcomes by focusing specifically on the part of the pathway mutated in UM. Activation of both downstream pathways, ERK and YAP, by oncogenic driver mutations in the CYSLTR2/GNAQ/GNA11/ PLCB4 pathway has been identified as critical for the clinical pathogenesis of UM, and combinations of inhibitors targeting these pathways are currently in early phase clinical trials50. One key aspect of PLCβ that we include in our model is that, in addition to acting as a Gαq/11 effector protein that can transmit signals downstream, PLCβ can also act as a GAP for Gαq/1151. Lastly, while the time dynamics of signaling pathways are critical in certain circumstances, here we focus on the steady-state activation of effectors because in cancer it is the sustained levels of growth signals that seem to be most relevant. Thus, we focus on the dynamic equilibrium that occurs at steady state. As we do not explicitly model signaling downstream of PLCβ and TRIO, we implicitly assume that there is a monotonic relationship between each of these effectors and their downstream pathways (i.e., increased TRIO activation results in increased Rho/Rac/FAK/YAP signaling, and increased PLCβ activation results in increased IP3/DAG/PKC/ERK signaling). Similarly, the model assumes a small basal level of receptor in the active conformation52 that implicitly takes the steady-state levels of β-arrestin desensitization into account. While the activation of a GPCR by its ligand can involve a complex set of intramolecular reconfigurations and is an important aspect of canonical G-protein signaling in response to time-varying ligand concentrations, here we focus on the steady-state activation levels relevant in oncogenic signaling. Thus, the main functional outputs for the model are the steady-state levels of activated TRIO and PLCβ (by either WT or mutant Gα subunits).
In summary, the model presented here attempts to capture the following key underlying biochemical steps: basal active CysLT2R catalyzes the exchange of GDP for GTP in the alpha subunits (Gαq/11) of heterotrimeric G-proteins downstream. Once bound to GTP, the Gαq/11 subunits effectively disassociate from the beta-gamma subunits (Gβγ). Once disassociated, the GTP-bound Gαq/11 can either directly bind and activate one the downstream effectors PLCβ or TRIO, disassociate from GTP to become nucleotide-free Gαq/11 or catalyze the hydrolysis of GTP to GDP through either intrinsic GTPase activity or through the catalyzed action of RGS. The cycle is then completed when GDP-bound Gαq/11 again binds Gβγ to form the heterotrimeric G-protein. The mutant CysLT2R and Gαq/11 proteins are modeled explicitly as distinct species in the model that undergo the same reactions described here but with different biochemical properties compared to the WT proteins based on available data.
The mathematical model reproduces key signaling outputs for different modeled disease genotypes
As a first evaluation of the mathematical model described above, we simulated the common Gαq/11 mutants: Q209L and Q209P. The Q209L/P mutants were considered heterozygous (i.e., 75% of Gα subunits are modeled as WT and 25% mutant, as described above). These assumptions are based on the TCGA data from cBioPortal which show that GNAQ Q209L/P mutations are typically heterozygous53,54. The model naturally and emergently results in significantly higher steady-state levels of complexes between activated (GTP-bound) Gα and its effectors TRIO and PLCβ (hereafter, “activated TRIO” and “activated PLCβ”) for each of the mutant Gαq cases when compared to the basal WT case (Fig. 2a). This suggests that the mechanisms considered by our mathematical model, when combined with approximate parameters for the WT and mutant proteins, are sufficient to explain the activation of these mutants.
a Active TRIO and PLCβ levels resulting from simulation of the mathematical model for cases where all GNAQ/11 is WT, where one GNAQ/11 allele is the oncogenic Q209L, and where one GNAQ/11 allele is the oncogenic Q209P allele for modeled basal conditions (no CysLT2R ligand-mediated activation) and for modeled stimulated conditions (CysLT2R ligand-mediated activation). b Active TRIO and PLCβ levels resulting from simulation of the model for modeled heterozygous (one copy) and homozygous (two copies) of GNAQ/11 Q209L, GNAQ/11 Q209P, GNAQ/11 R183C, and CYSLTR2 L129Q compared to the modeled all GNAQ/11 and CYSLTR2 WT genotype.
We considered the anticipated behavior if the oncogenic mutation was homozygous instead of heterozygous. Modeling this case, where there was now twice as much of the mutant protein abundance in the model but the same total amount of wild-type plus mutant Gαq/11, resulted in a slight increase in total steady-state signal of activated TRIO and activated PLCβ (Fig. 2b). We extended our model to include the less strongly activating R183C Gαq/11 mutations. Introduction of the approximate mutant parameters also resulted in constitutive signaling of the pathway for this mutant in both heterozygous and homozygous cases (Fig. 2b).
The mathematical model reproduces unanticipated responses to pharmacological G-protein inhibition
One of the major goals of understanding the molecular causes of tumorigenesis in UM is to develop better therapies. Most efforts in the clinic up to this point have focused on direct inhibition of downstream effectors such as MEK and PKC or on an immunotherapy approach. Both general approaches have had limited success39.
One key development in the field over the last several years has been the discovery and exploration of direct pharmacological inhibitors of Gαq/11 subunits such as YM-254890 (YM) and the naturally occurring cyclic depsipeptide FR900359 (FR)35,36,37,38. Specifically, FR has been shown to inhibit Gαq/11 signaling in UM cells and mouse xenograft models37,38. The current consensus is that FR acts as a guanine nucleotide disassociation inhibitor (GDI), effectively locking Gα q/11 in the GDP-bound, inactive state, although it is unclear if this is the only mechanism responsible for the inhibition of oncogenic signaling observed in various in vitro and in vivo models of UM38. While it has been recently shown that the major effects of FR/YM are in fact on-target for Gαq/1155, the efficacy of a drug that targets the inactive (GDP-bound) form of a mutant protein that is “constitutively active”, where the active form is GTP-bound, has created confusion about this mechanism of action. This confusion regarding mechanism is similar to the issues originally surrounding the direct targeting the KRAS G12C oncoprotein with inhibitors that target the GDP-bound form of the constitutively active KRAS G12C mutant9,56,57,58.
As another test of whether mathematical modeling can be useful in situations that are confusing to domain experts of the biological subject matter, we expanded the model to include drug targeting. Specifically, we modeled a FR/YM-type GDI drug that can reversibly bind and sequester GDP-bound Gα subunits (Fig. 3a). We then simulated the model at varying drug concentrations for both WT/WT and WT/mutant cases (Fig. 3b). The model suggests the proposed FR mechanism of binding to the GDP-bound form of the “constitutively active” Gα mutants is sufficient to explain its action. In other words, modeling reveals that the issues that superficially appeared confusing are not actually problematic and do not indicate a gap in mechanistic understanding. Interestingly, the model is able to replicate a reported59 log-shift in sensitivity between Gαq mutants that are known to be partially GAP-sensitive (e.g., R183Q/C)60 and more common oncogenic mutants such as Gαq Q209L (Fig. 3b). The model’s ability to predict relative differences in response to pharmacological inhibition, including the ability to reproduce the non-intuitive effects of a drug that binds selectively to GDP-bound forms of Gαq/11, help validate our mathematical model and suggest it may have future value to help contextualize and interpret drug development studies for UM.
a Schematic of extensions to the model presented in Fig. 1b to include a Gαq/11-inhibitor. The schematized G-protein alpha subunits (red circles) can represent either (or both) wild-type and mutant forms. b Model simulations of Gαq/11-inhibitor dose responses on modeled GNAQ/11 WT, heterozygous R183C, and heterozygous Q209L genotypes.
Mathematical model illuminates unknown differences in downstream signal activation between mutant Gαq and mutant CysLT2R
UM cases without GNAQ/GNA11 mutations often contain activating mutations upstream in the associated receptor CysLT2R61,62 or in the canonical downstream effector of Gαq, PLCβ4 (the protein encoded by PLCB4)63 (Fig. 4a). Most of the genomic data on UM is for GNAQ/GNA11 mutant cases and most clinical data is from patients with a GNAQ/GNA11 mutation. Additionally, while there are several preclinical experimental model systems to study GNAQ/GNA11 mutant UM, from cell lines to mouse models64, there are no equivalent experimental model systems for CYSLTR2 mutant UM65. Although the biochemical properties have been measured for the most common point mutations, Q209L and Q209P, in Gαq/1144, similar biochemical data for point mutations in CysLT2R (the upstream receptor that activates Gαq) are not available. We hypothesized that, despite the gaps in knowledge, our mechanistic model mapping between biochemical properties and cellular phenotypes may be able to clarify ranges of potential behaviors and be able to generate insights that prioritize new areas for experimental and/or bioinformatic investigation.
a Pie chart of mutations in the pathway involving CYSLTR2/GNAQ/GNA11/PLCB4 for 80 UM patients reported in previous studies53,54. b Schematic of oncogenic Gαq and oncogenic CysLT2R signaling in the model. The schematized G-protein alpha subunits (red circles) can represent either (or both) wild-type and mutant forms. c Sobol total sensitivity index for the difference in active TRIO and PLCβ between the GNAQ Q209L and CYSLTR2 L129Q mutant settings based on simulation of the mathematical model. See the supplement for details on the parameters varied and ranges. d Schematic highlighting two diverging behaviors suggested from the extended simulations of the model across parameter space: (1) that CysLT2R mutant networks are less able to activate TRIO/YAP/TAZ than Gαq/11 mutant networks, or (2) that CysLT2R mutant networks are less able to activate PLCβ/PKC/ERK than Gαq/11 mutant networks. Which situation occurs will depend on where the actual system is within parameter space, and could be determined by experimentally measuring parameters and/or by experimentally measuring cellular phenotypes. The bar charts on the right display representative sample cases of the two diverging behaviors; fold change in activation relative to the WT case is used as a measure of pathway activation within these bar charts.
We updated our model to the most common CysLT2R mutation: L129Q (Fig. 4b). This mutation was taken to be heterozygous based on TCGA data and modeled as a higher surface concentration of activated receptor (i.e., 50% of the CysLT2R receptor population was modeled as WT and 50% mutant). Recent work has demonstrated that the CysLT2R L129Q mutant poorly recruits β-arrestin to avoid desensitization66. Current opinion in the field is that the mutually exclusive GNAQ/11 and CYLSTR2 mutations provide similar oncogenic signals, and our initial parameterization estimates based on the best available data for the CysLT2R L129Q mutant resulted in an elevation of active TRIO that was slightly to moderately weaker than what the model predicts for Gαq/11 mutations (Fig. 2b). However, we also noticed that our initial parameter estimates resulted in a near complete loss of PLCβ activation compared to the Gαq/11 mutations (Fig. 2b). We hypothesized that this might suggest our initial parameter estimates need to be adjusted.
To further explore the potential ranges of behaviors with the model for alternative parameter sets, we evaluated the potential similarity of model-predicted behaviors for CysLT2R and Gαq/11 mutant conditions over a wide range of concentrations of the proteins considered in the model and over a wide range of kinetic rate constants for the considered reactions. For each parameter set, we simulated our mathematical model of the Gαq/11 Q209L and CysLT2R L129Q mutations. We then performed a global sensitivity analysis to compare the downstream activation of these two mutations over the whole parameter range (Fig. 4c and the supplement). The rate constant for GAP-stimulated hydrolysis of Gαq/11-bound GTP by PLCβ (denoted by khyd2) was determined to be the most sensitive parameter with respect to the difference in PLCβ activation between the two mutants, while the total concentration of TRIO was the most sensitive parameter with respect to TRIO activation. As the biochemical parameters of this system have not been fully measured, this sensitivity analysis is evaluating the scope of potential behaviors that are within the realm of possibility of available knowledge as much as it is evaluating the robustness of the model.
Although the prevailing assumption is that GNAQ/11 and CYSLTR2 mutant networks are effectively equivalent, we were intrigued that our model suggested that the available data do not require them to be equivalent (Fig. 4d). The unknown parameters of the system could potentially allow CYSLTR2 mutant UM to have strong ERK but weak YAP signaling (~5% of parameter sets sampled), to have weak ERK but strong YAP signaling ( ~ 50% or parameter sets sampled), or to have weak activation of both pathways (~25% of parameter sets sampled). While it is difficult to translate the relative frequency of behaviors over the entire parameter range into predictions about the true parameter values and behavior, we were intrigued at the possibility of there being a lack of functional equivalency between the GNAQ/11 and CYSLTR2 mutants as this is counterintuitive based on the canonical view of the pathway. As recent studies have highlighted both ERK and YAP signaling as being critical to GNAQ/11 mutant UM, we thought it was important to experimentally determine whether CysLT2R L129Q can also activate both pathways.
Experiments reveal that while oncogenic GNAQ can activate ERK and YAP signaling, oncogenic CYSLTR2 can only strongly activate ERK
To experimentally explore model-suggested potential for divergent cellular phenotypes, we expressed WT or mutant GNAQ or CYSLTR2 constructs in HEK 293 T cells and quantified downstream signaling and activation of ERK and YAP signaling via RNA-sequencing. HEK 293 T cells have often been used to explore signaling from these mutants48,49,61,67 and provide a convenient system for comparing signaling phenotypes of the mutants. (Note: UM cell lines that allow comparisons between GNAQ/GNA11 and CYSLTR2 mutants in an isogenic background are not currently available). Gene-set enrichment analysis (GSEA)68,69 revealed that both mutants yielded a statistically significant enrichment in the hallmark “KRAS signaling up” signature (active KRAS activates ERK) and the “YAP conserved” signature70 compared to the mock transfection control (Fig. 5a, Supplementary Fig. 1). In addition, when the GNAQ Q209L expressing cells’ gene expression signature was compared to the gene expression signature of the cells expressing CYSLTR2 L129Q, it was noted to be significantly enriched in the “YAP conserved” signature (suggesting that there was much more YAP signaling in the GNAQ Q209L cells) while there was no enrichment in the hallmark “KRAS signaling up” signature (suggesting that both GNAQ and CYSLTR2 mutants activated KRAS/ERK to more similar levels) (Fig. 5b). We validated these signaling differences by western blot (Fig. 5c). Based on these experiments, the CYSLTR2 L129Q mutant appears to have strong ERK activation but weak YAP signaling, in contrast to the GNAQ Q209L mutant, which seems to activate both pathways strongly (Fig. 5d).
a Gene set enrichment analysis (GSEA) of differentially expressed genes for GNAQ Q209L and CYSLTR2 L129Q transfected 293T cells, each compared with mock transfected cells. b GSEA of differentially expressed genes between the GNAQ Q209L and CYSLTR2 L129Q transfected conditions for the same cells as in (a). c Total cell lysates were probed by western blot for phospho-ERK, phospho-YAP, phospho-FAK, Gαq, and GAPDH. Blots are representative of n = 3 experimental replicates. d Schematic of the experimentally observed differences in the relative strength of downstream activation for the Gαq Q209L and CysLT2R L129Q mutations, which identifies which of the two possible diverging behaviors described in Fig. 4d is experimentally demonstrated to occur here. The schematized G-protein alpha subunits (red circles) can represent either (or both) wild-type and mutant forms.
Iterative evaluation and adjustment of model parameters based on the new experimental data
We determined which simulated alternative parameter sets matched our experimentally observed activation patterns and which did not (i.e., which parameter sets gave the result that the CysLT2R L129Q mutation yielded greater than or equal levels of PLCβ activation and lower levels of TRIO activation compared to the Gαq Q209L mutation). The two parameters which showed the most significant discrimination between which simulations matched the experiments and which simulations did not match the experiments based on a Kolmogorov-Smirnov test were: the bias of the Gαq Q209L mutant binding effectors to PLCβ relative to TRIO (denoted by the ratio of the association rate constants for each reaction: kat2mut/wt/kat1mut/wt) and the rate constant for GAP-stimulated hydrolysis of Gαq-bound GTP by PLCβ (khyd2) (Fig. 6a, Supplementary Fig. 2 and the supplement). A receiver operating characteristic (ROC) analysis suggested that, for the given parameter ranges, kat2mut/wt/kat1mut/wt was the more effective binary classifier for which simulations matched the experiments and which simulations did not match the experiments (Fig. 6b). Of note, co-immunoprecipitation studies that compared the abilities of Gαq Q209L and Gαq WT to bind to TRIO and PLCβ suggested a modest bias for TRIO over PLCβ effector binding for Gαq Q209L compared to Gαq WT44. The rate constant for GAP-stimulated hydrolysis of Gαq-bound GTP by PLCβ (khyd2) was another critical parameter. Parameter sets in which khyd2 was ~1–40-fold faster than basal hydrolysis were more likely to match the experimental results than parameter sets in which khyd2 was >50-fold faster than basal hydrolysis (Fig. 6a).
a Normalized distributions of parameter values: the bias of the Gαq Q209L mutant binding effectors (denoted by kat2mut/kat1mut) and the GAP-stimulated hydrolysis of Gαq-bound GTP by PLCβ (denoted here as the fold over basal hydrolysis khyd2/khyd) for simulations that qualitatively reproduced the experimental activation patterns (shown in pink) and simulations that did not (shown in cyan). A total of n = 11,200 parameter sets were simulated (Fig. 4c and the supplement). b Plotted are the receiver operating characteristic (ROC) curves for using either: (1) the bias of the Gαq Q209L mutant binding effectors (denoted by kat2mut/wt/kat1mu/wt) with a threshold between kat2mut/wt/kat1mu/wt ≥1 and kat2mut/wt/kat1mu/wt ≥14, or (2) the rate constant for GAP-stimulated hydrolysis of Gαq-bound GTP by PLCβ (denoted here as the fold over basal hydrolysis khyd2/khyd) with a threshold between khyd2/khyd <1 and khyd2/khyd <100, as a binary classifier to identify whether the model simulated at various parameter sets either matched or did not match the experiment. It is worth noting that khyd2/khyd was set to roughly 770 in Fig. 2 based on the available data. c Active TRIO and PLCβ levels resulting from simulation of the mathematical model for the WT, GNAQ Q209L, and CYSLTR2 L129Q settings over a range of GAP-stimulated hydrolysis of Gαq-bound GTP by PLCβ (denoted here as the fold over basal hydrolysis khyd2/khyd). The pink region denotes the parameter range for which the simulations qualitatively reproduced the experimental activation patterns and the cyan region denotes the parameter range for which the simulations did not qualitatively reproduce the experimental activation patterns. Based on the experimental observations in Fig. 5 and the analysis in (a, b) above, the bias of Q209L TRIO/PLCB binding was set to kat2mut/wt/kat1mu/wt = 4 for these simulations.
To explicitly compare the behavior of the Gαq Q209L and CysLT2R L129Q mutations at varying values of khyd2, we plotted the levels of active PLCβ and active TRIO for each (Fig. 6c), while keeping other parameters fixed. This suggests that when PLCβ has very strong GAP activity (i.e., rapidly promotes the conversion of GTP-bound Gα subunits to GDP-bound Gα subunits), our model of the CysLT2R L129Q mutation will not yield significant steady-state activation of PLCβ because all of the Gα subunits are WT and thus fully GAP sensitive (i.e., any Gαq-bound GTP will be quickly hydrolyzed by PLCβ resulting in termination of the signal). Overall, these experimental and computational analyses suggest that there is a significant difference in the downstream activation patterns of activating GNAQ and CYSLTR2 mutations. Specifically, the CYSLTR2 L129Q mutant appears to be relatively deficient in activating the TRIO- > FAK, YAP pathway compared to the GNAQ Q209L mutant.
We also present an updated set of parameter estimates that is similar to our initial parameters, but where we have chosen new values for parameters kat2mut/wt/kat1mut/wt and khyd2 based on this analysis (see the supplement for details). These new parameters are able to reproduce all of the same behaviors as our initial parameter set and they also better reproduce the qualitative differences in downstream signaling observed during our new experiments with CYSLTR2 and GNAQ mutants (Supplementary Fig. 3). We therefore suggest these updated parameters, combined with the model, may serve as a good foundation for further model-driven studies of UM.
Bioinformatic analysis of patient data reveals CYSLTR2 mutations co-occur with semaphorin/plexin gene mutations
Based on our computational and experimental analysis presented in the previous sections, we propose the hypothesis that CysLT2R L129Q mutations may be deficient in TRIO- > FAK, YAP pathway activation compared to Gαq Q209L mutations. This hypothesis raises the important question: what is happening in UM patients with CYSLTR2 L129Q mutations? While there are many potential explanations, one intriguing possibility is that there are compensating mutations in UM patients with CYSLTR2 L129Q mutations that provide additional FAK- > YAP pathway activation (Fig. 7a). None of the UM genomics studies have reported such co-occurring mutations, but the limited number of CYSLTR2 mutant UM (less than ten cases sequenced to date) leaves this subset underpowered for unbiased discovery analysis of the mutation data. Of note, in a previous bioinformatics analysis of TCGA data, CYSLTR2 mutant UM patients do not localize to specific molecular or clinical subsets of UM43, suggesting there is ultimately similar signaling activity in these patient tumors. We also performed a differential expression analysis of GNAQ/GNA11 mutant vs. CYSLTR2 mutant patients in TCGA followed by GSEA, which revealed no statistically significant difference in either of the “KRAS signaling up” or “YAP conserved” signatures mentioned previously (Supplementary Fig. 4a). It is worth noting that of the 14 genes identified as significantly differentially expressed between GNAQ/GNA11 and CYSLTR2 mutant UM patient samples in TCGA, PRAME (preferentially expressed in melanoma) has previously been identified as independent marker of metastasis in UM71 (Supplementary Fig. 4b).
a With CysLT2R L129Q deficient at activating FAK/YAP/TAZ signaling relative to Gαq Q209L, we hypothesize that there may be co-occurring genetic and/or epigenetic changes that contribute to FAK/YAP/TAZ activation. The schematized G-protein alpha subunits (red circles) represent the wild-type form. Gene Ontology (GO) analysis of the mutations observed in the four published, exome sequenced, CYSLTR2 L129Q UM53,54 identifies the Axon Guidance Pathway as the most enriched pathway. b Schematic highlighting 2 of the 3 patients with CYSLTR2 L129Q mutation in the TCGA have a semaphorin/plexin co-mutation, which was done in our discovery phase. We then investigated a separate Uveal Melanoma data set (ref. 82) and found two of the two CYSLTR2 L129Q patients had a semaphorin or plexin gene mutated. c Schematic of assay to treat UM 92.1 cells with the selective FAK inhibitor VS-4718 and collect samples for bulk RNA-sequencing. d Percentage of genes that were significantly up or downregulated on all 4 days of sample collection genome-wide vs. the semaphorin/plexin pathway. e Gene Ontology (GO) analysis of the 36 genes that were significantly up or downregulated on all 4 days of sample collection.
Although the number of genomically sequenced CYSLTR2 mutant UM samples may be underpowered for unbiased discovery analyses, we speculated that the number may be sufficient for focused exploration of genes that could potentially activate FAK- > YAP signaling. We first analyzed UM patient data in TCGA using cBioPortal53,54. A gene ontology (GO) analysis72,73,74 of the list of genes mutated in these CYSLTR2 mutant patients using the Reactome Knowledgebase75 yielded a statistically significant enrichment of the axon guidance pathway (Fig. 7a, Supplementary Fig. 4c). Inspection of the genes revealed a statistically significant enrichment of mutations in the semaphorin/plexin signaling family in UM patients with CYSLTR2 L129Q mutations (Fig. 7b). Semaphorins are membrane-bound or diffusible ligands that signal through the plexin receptors and were originally identified as important in axon guidance and angiogenesis but have more recently been implicated in several aspects of cell-cell communication, cancer76 and FAK activation77,78,79. Previous work has identified that BAP1 loss in UM can lead to significant deregulation of axon guidance pathways, including several semaphorin/plexin genes80 and another study identified several semaphorin genes as differentially expressed in melanocytes expressing GNAQ Q209L81. Thus, there were some reasons to speculate that the observed mutations in the semaphorin/plexin pathway may be providing the hypothesized contribution to FAK activation. Importantly, UM is relatively unique among malignancies that impact adults in that there are typically a very low number of coding mutations. Thus, there were few co-occurring mutations within the CYSLTR2 mutant cancers (between 5 and 19 coding mutations per sample), intuitively reducing the chance that the enrichment of semaphorin/plexin mutations in CYSLTR2 mutant UM relative to GNAQ/GNA11 mutant UM was a false discovery. Our analysis of co-occurring BAP1 monosomy or mutation as a covariate found no significant correlation with CYSLTR2 mutations and BAP1 monosomy/mutation.
To prospectively test whether this observed enrichment of mutations existed in other available datasets, we found and analyzed additional UM genomics data82. Of note, these data included two additional patients with CYSLTR2 L129Q mutations and both samples have coding mutations in the semaphorin/plexin pathway (Fig. 7b). Thus, a total of 4/5 patients with CYSLTR2 L129Q mutations from TCGA and ref. 82 have at least one co-occurring mutation in the semaphorin/plexin signaling family. Analysis of the observed semaphorin/plexin mutations observed in TCGA and ref. 82 with PolyPhen283 revealed two of the mutations (PLXNA4 T642I and PLXNA4 R1626Q) were predicted to be damaging while two mutations (PLXND1 Q657H and SEMA7A E183K) were predicted to be benign. The observed mutations in PLXNA4 and PLXND1 in TCGA are found in similar domains of the respective proteins (Supplementary Fig. 4d). We also observed a statistically significant difference in disease-specific survival in UM patients with low expression of PLXNB1 and PLXNA1 mRNA and high expression of PLXNC1 and SEMA4D mRNA, thus providing additional circumstantial evidence that plexin/semaphorin signaling has a role in uveal melanoma (Supplementary Fig. 5a). There are several reports in the literature of the tumor suppressor effects of semaphorin/plexin signaling in cutaneous melanoma84,85,86. Thus, although a role for semaphorin/plexin signaling in uveal melanomagenesis has only been briefly mentioned in the literature80,81, our analysis finds a variety of suggestive data.
Treatment of GNAQ mutant uveal melanoma cells with FAK inhibitors further suggests a role for semaphorin/plexin signaling in UM
FAK inhibitors are being evaluated for use in UM50,87. We hypothesized that the semaphorin/plexin pathway may play a role in the response to FAKi targeted therapy in the UM context. To experimentally test this hypothesis, we continually passaged GNAQ mutant UM 92.1 cells in either 100 nM of the selective FAK inhibitor VS-4718 or DMSO control for 21 days (Fig. 7c). At this dose, there is a strong suppression of signaling (Supplementary Fig. 5b) but little to no suppression of proliferation (Supplementary Fig. 5c). We hypothesized that there may be compensatory changes in semaphorin and/or plexin genes that counteract the targeting of FAK to maintain proliferation at these modest doses of FAK inhibitor. To evaluate, we performed RNA-sequencing on samples collected at 0, 3, 7, 14, and 21 days. Differential expression analysis revealed that several of the semaphorin/plexin genes were significantly up or downregulated in response to continued FAK inhibitor treatment compared to DMSO control, with 3/29 or 10.34% of these genes significantly upregulated on all 4 days of sample collection (Fig. 7d). In comparison, the background rate of genes significantly up or downregulated on all 4 days of sample collection was 0.21%. A GO analysis of the 36 genes that were significantly up or downregulated on all 4 days of sample collection revealed that the semaphorin-plexin signaling pathway was the most significantly enriched GO term and that all of the GO terms that were significantly enriched contained SEMA6A and SEMA6D (Fig. 7e). The patterns of expression for each semaphorin/plexin gene in response to FAK inhibitor treatment can be seen in Supplementary Fig. 5d.
Discussion
In this study, we investigated whether biochemical-mechanism based mathematical models of critical oncogenic signaling networks could be used to help open new areas for research into mechanisms of disease.
We here found this to be the case: the model-based analysis first revealed that pathogenic and mutually exclusive GNAQ/11 and CYSLTR2 mutations are not functionally equivalent. We found that although the dual-activation of ERK and YAP pathways that is essential for UM can be generated by a canonical GNAQ mutation, a pathogenic CYSTLR2 mutation is only able to strongly activate the ERK pathway and may require an additional, co-occurring event to activate the YAP pathway. The lack of equivalency was first suggested by our mathematical modeling, tested and confirmed experimentally, and then evaluated with genomic and transcriptomic analyses. This new understanding uncovered the potential mechanism that CYSLTR2 mutant UM use to activate the YAP pathway: deregulation of the semaphorin/plexin pathway.
There are reports of various semaphorin/plexin mutations in the context of cutaneous melanoma and other cancers88,89. However, a role for semaphorin and plexin signaling in UM does not have extensive prior evidence. Our mathematical modeling led us to uncover signaling differences between oncogenic Gαq and oncogenic CysLT2R, which we then experimentally confirmed. We more generally hypothesize that there may be signaling differences for PLCB4 mutant melanoma and for other mutant forms of GNAQ/11 and of CYSLTR2. Additional characterization of signaling by other mutants found in UM appears to be an important area for more experimental investigation that our model-based study has motivated.
Our hypothesis that secondary, compensating mutations are required for CYSLTR2 mutant UM provides a potential explanation for the relative scarcity of UM patients with CYSLTR2 L129Q mutations. If both Gαq/11 and CysLT2R mutations behaved exactly the same in terms of signaling activation, they might be expected to occur in more even proportion in UM patients. On the other hand, if CYSLTR2 L129Q mutant tumors need a second compensating mutation to become fully oncogenic, it is logical that they would be rarer. An important test of our hypothesis in the future will be if the enrichment of semaphorin/plexin mutations we identified for CYSLTR2 mutant UM patients is confirmed as more patient samples are sequenced. Of note, there was one additional UM patient sample90 that harbors a CYSLTR2 L129Q mutation and that is not in either of the datasets mentioned previously. This patient was identified to have a PLXND1 coding mutation90, but upon our re-analysis of the raw data we determined that the mutation was intronic. While we cannot easily determine the meaning of this mutation, it is at the very least an intriguing coincidence, and could potentially be functional91,92. Another area of future exploration suggested by this model-driven study is additional sequencing of UM patients, specifically to evaluate co-occurrence of plexin and semaphorin gene mutations. Biochemical and biophysical characterization of plexin and semaphorin gene mutations is another emerging area for future study, as motivated by this study.
The enrichment of plexin and semaphorin mutations in CYSLTR2 mutant UM patients and evidence of survival differences based on gene expression make a strong case that plexin and semaphorin proteins play a role in uveal melanoma. As this is a large family with related, but diverging, roles, a more thorough characterization of the different plexins and semaphorin proteins in uveal melanoma, including context-specific function, cross-talk, and feedback would provide important information that is needed to more fully interpret experimental data and could potentially identify new therapeutic strategies.
Our computational analyses reveal that the effector binding properties of the Gαq mutant and the rate constant for GAP-stimulated hydrolysis of Gαq-bound GTP by PLCβ are critical for our mathematical model to reproduce the experimentally observed activation patterns. This result was not intuitive to experimental biologists working in the field or to computational biologists working on the model. That the behavior of the system depends upon the specific parameters of the model further highlights that it is not possible to simply look at a model schematic and infer model outcomes. However, once an empirical observation is made and parameters of the model are found that reproduce the empirical observations, the model may facilitate the development of an intuitive understanding of the system behavior. In this study, we experimentally found that CysLT2R L129Q is impaired at activating TRIO/YAP signaling relative to Gαq/11 and does not seem to be impaired at activating PLCβ/ERK. This was the opposite of the computational inference that resulted from using our initial parameter estimates.
Our revised parameter estimates are consistent with the experimental observations and consideration of how parameters needed to change to match the observations does help produce an intuitive understanding of the system. The observed behavior of the two mutants, and the parameter changes, may be intuitively understood as follows: (a) our initial estimate for PLCβ-mediated hydrolysis of active, GTP-bound, WT Gαq/11 was too large, resulting in a profound inhibition of PLCβ (and thereby ERK) activation even when CysLT2R L129Q was driving WT Gαq/11 activation. Correspondingly, alternative parameter sets where the value of PLCβ-mediated hydrolysis of active, GTP-bound, WT Gαq/11 was smaller than the original estimate were more likely to result in activation of PLCβ/ERK. The other parameters that needed to be updated for the model to be more consistent with empirical observations involved the association rate constants between mutant Gαq/11 with its two effectors, TRIO and PLCβ. If Gαq/11 Q209L has similar affinity for TRIO as WT Gαq/11, and if both Gαq/11 Q209L and CysLT2R L129Q induce similar total levels of GTP-bound Gαq/11 (i.e., WT and mutant combined), then one would expect Gαq/11 Q209L and CysLT2R L129Q to bind TRIO and activate YAP/TAZ signaling similarly. However, the experiments found that CysLT2R L129Q was impaired at activating YAP/TAZ signaling. Our updated parameters therefore assume an elevated association rate constant for Gαq/11 Q209L binding to TRIO so that the model better matches the experimental observation.
Our finding that CysLT2R L129Q is not impaired at activating ERK suggests that the PLCβ isoforms present in uveal melanoma catalyze hydrolysis to a more modest extent. A consideration of experimental reports of these parameters can help evaluate the potential validity of this inference. For our initial parameterization, we chose the rate constant for GAP mediated hydrolysis by PLCβ based on one study that reported the rate constant for GAP mediated hydrolysis of PLCβ was up to 1000-fold faster than basal hydrolysis for PLCβ151. Other isoforms in other studies are reported to have a wide range of rate constants93. The rate constant for GAP mediated hydrolysis of PLCβ in our model can be thought of as an aggregate effect from all the PLCβ isoforms present. It should also be noted that the binding of Gαq to PLCβ has been reported to exhibit anomalous affinity94 and there may be more complex activation mechanisms that are beyond the scope of our current model.
We highlight several important qualifiers for our parameter sets. Our initial parameter estimates represented reasonable first estimates based on the literature, such as what one would use when initiating a new study. Our revised parameter estimates are an updated set that is more consistent with the experimental observations, but is not formally fit or optimized. If and when parameters are formally fit to a model, it is common for there to be many alternative parameter sets that can reproduce the behavior of a system95. We could therefore have potentially identified a completely different alternative parameter set that would have been effectively equally consistent with the behaviors we would like the model to match. Therefore, our updated parameters should not be interpreted as improved estimates for the underlying biochemical parameter. As another way to state this, in our computational analysis we kept the equilibrium dissociation constants for protein-protein interactions constant and only allowed the association rate constant to vary. Had we instead varied the dissociation rate constant while keeping the equilibrium dissociation constant fixed, we would have estimated changes to dissociation rates that would yield the same overall change in equilibrium dissociation constant. Overall, many other alternative parameterizations could have been found. Future studies could attempt to better constrain parameters, and these overall results highlight a thorough biochemical characterization of PLCβ as an important area for future research in the study of UM molecular pathogenesis. Additionally, the empirical observation of CysLT2R L129Q is not impaired at activating ERK but is relatively impaired at FAK activation further suggests that Gαq/11 has a bias for binding to PLCβ relative to TRIO; this is another model-based inference that can be pursued in subsequent experimental work.
Lastly, it should also be noted that there are several possible future extensions to the model that could complement any future experimental studies. As an example, explicitly including components of the semaphorin/plexin pathway in the model could help explore possible mechanisms of co-activation with the CysLT2R L129Q mutant receptor. Although several aspects of semaphorin and plexin signaling have been investigated with computational and mathematical models96,97,98,99,100 biochemical-mechanism based mathematical models do not appear to have heavily utilized with respect to investigations of signaling from plexins through its downstream effectors. Mathematical models that investigate cross-talk between Gαq/11 and semaphorin/plexin signaling could potentially be enlightening; for example, a simple model may help estimate how much plexin-driven FAK activation would be needed in combination with a CysLT2R L129Q mutant to achieve a total signal comparable to what a Gαq Q209L would induce. We envision this as an area where mathematical modeling could be particularly useful. Other potential modeling directions include explicitly modeling downstream from PLCβ and TRIO, as well as expanding the characterization of PLCβ to investigate pathogenic PLCβ mutations.
One challenge in the development of a mechanistic model is the identification of the appropriate parameter values. For example, even though GPCR and G-protein signaling has been extensively studied, the characterization of fundamental rate constants for reactions and for mutants is an ongoing area of research44,101,102,103,104. Thus, published mechanistic models that focus on GPCR signaling and downstream G-proteins must estimate and approximate parameter values. In contrast, the RAS GTPases (HRAS/NRAS/KRAS) have been very well studied at the level of component reactions, including for several oncogenic mutations14,15,105. This has enabled the development of mathematical models of RAS signaling that can directly incorporate these experimental data as parameters7,106,107. These models of RAS have been useful for understanding why pathogenic mutations are activating, for inferring new aspects of cancer genetics, and for understanding mutation-specific treatment responses. However, it has been unclear whether biochemically detailed mechanistic models can also add a meaningful contribution to studies of pathogenic mutations when the chemical-kinetic properties of the pathway and mutations are far less thoroughly quantified. Our work here demonstrates that this approach can still catalyze progress and uncover previously unknown features of a system. It is worth noting that mutations in the CYSLTR2/GNAQ/GNA11/PLCB4 signaling pathway have been identified in other malignancies such as blue nevi108,109 and tumors of the central nervous system110,111,112,113, and somatic GNAQ mutations are known to cause the congenital disorder Sturge-Weber syndrome67. Looking forward, our mathematical model could be applicable to these cases as well.
This study highlights how model-based inference can help break new ground through the re-analysis of existing data to suggest new ideas that can be tested with focused, hypothesis-driven, experiments. We refer to this approach as MAGPIE (Model Assisted Generation of Predictions and Interpretation of Experimental data) (Fig. 8). The value of this approach can be better understood if one considers why some of these observations were not previously made. For example, although our analysis relies upon readily available genomic data, previous analyses of these data did not identify the co-occurrence of semaphorin/plexin pathway mutations in CYSLTR2 mutant UM43,114. This is likely because a relatively small number of UM patient samples have received exome-level sequencing; only 106 cases are currently listed in cBioPortal. Of those, 97 have GNAQ/GNA11 mutations and nine do not. This is an extremely small dataset, and a non-focused analysis that considers all pathways is unlikely to yield meaningful results. Our mechanistic mathematical modeling led to experimental discoveries, which in turn suggested very focused analyses of the limited genomic data. In this manner, we were able to uncover a new aspect of UM biology that was featured in the available data (biochemical, biophysical, and genomic) but that may not have been detected without an integrated, mechanistic approach to evaluate the pathway. It should be noted that there are several prominent examples in the literature utilizing the combination of mechanistic modeling and data-driven approaches to generate insight into diverse biological systems115,116.
This work demonstrates how an integrated approach involving mechanistic modeling, bioinformatic analysis, and experimental approaches can uncover new directions and aid in our understanding of cancer biology. The figure was created by illustrator Amy Cao with express permission to use.
Our work highlights the power of integrating multiple computational and bioinformatics methods with experimental and genomic biology. Such integrated studies that iteratively alternate between (new or existing) experimental data and mathematical and computational analyses have great potential for elucidating unknown aspects of cancer biology. For example, the integration of additional datasets to generate novel, evidence-based, hypotheses can allow for focused queries of genomic data where the signal of a biologically important process can be better uncovered than through an unbiased, all-hypotheses considered query of genomic data. Our ability to extract new insights into UM, a rare cancer with a small number of cases that have been sequenced, by focusing on an even smaller number of cases attests to the power of these integrated approaches as a tool for moving forward when vast datasets are not available and mechanism agnostic big data approaches cannot be utilized.
Methods
Mathematical model
A description of the mechanistic mathematical model of oncogenic Gα signaling presented here can be found in the main text and the explicit reactions, rate constants, protein abundances, and relevant references can be found in the supplement. Our model is a set of ODEs that we solve numerically using Python (RRID:SCR_008394). The global sensitivity analysis presented in Fig. 4c of the main text was performed using the SALib Python library117. Code that allows all computational results and figures to be reproduced is provided in the form of a Python notebook at https://github.com/StitesLab/GNAQ_model. A PDF that shows the executed Python notebook with outputs is also included for convenience.
Cell line models and culture method
HEK 293T cells were grown in DMEM supplemented with fetal bovine serum (FBS) (10%), penicillin (100 U/ml), and streptomycin (100 μg/ml). 92.1 UM cells (RRID:CVCL_8607) were grown in RPMI media containing 10% FBS (Atlanta Biologicals) and 1% penicillin/streptomycin (Corning).
Expression plasmid transfection
HEK 293T cells were plated in a 6-well plate in DMEM supplemented with 10% FBS, penicillin (100 U/ml), and streptomycin (100 μg/ml) 24 h before transfection. The following day, cells were transfected with expression plasmids with duplex containing 0.25 μg of DNA, 250 μl of optimem, and 10 μl of Lipofectamine 2000 per well. All constructs were provided by the Gutkind lab and mutations were verified by sequencing in the Stites lab. Cells were harvested 24 h after transfection.
Western blot analysis
Cell lysates were generated using Lysis buffer (Thermo Fisher Scientific, 1862301) containing protease inhibitor cocktail (Cell Signaling Technology) and incubated on ice for 1 h, with brief vortexing every 5 min. The total protein concentration was determined by Pierce Protein assay (Thermo Fisher Scientific). Protein samples were resolved by electrophoresis on 12% SDS–polyacrylamide gels and electrophoretically transferred to polyvinylidene difluoride (PVDF) membranes (Millipore Corporation) for 20 min at 25 V with the trans-blot turbo (Bio-Rad Laboratories). The blots were probed with the appropriate primary antibody and the appropriate fluorophore-conjugated secondary antibody. The protein bands were visualized using the Licor CLx Odyssey imaging station (Licor Biosystems) (RRID:SCR_014579). Comparative changes were measured with Licor Image Studio software from independent experiments. The antibodies used are: anti-phospho-Thr202/Tyr204-ERK1/2 (BioLegend, catalog #675502) (RRID:AB_2565604) (1:1000 dilution); anti-phospho-Ser127-YAP (Cell Signaling Technology, catalog #4911) (RRID:AB_2218913) (1:1000 dilution); anti-phospho-Tyr397-FAK (Cell Signaling Technology, catalog #8556) (RRID:AB_10891442) (1:1000 dilution);anti-GNAQ (Cell Signaling Technology, catalog #14373) (RRID:AB_2665457) (1:1000 dilution); anti-GAPDH (Santa Cruz Biotechnology, catalog #47724) (RRID:AB_627678). Multiple blots were run in parallel from the same lysates. Uncropped blots are provided as a source data file within the supplementary methods.
Bioinformatic analysis of somatic mutations in patient samples
To calculate the probability of the observed enrichment of semaphorin/plexin family mutations in TCGA (RRID:SCR_003193) and ref. 82 cited in the main text occurring by chance, we explicitly calculated the number of ways of having the observed number of co-mutations divided by the total number of ways of distributing the mutations at random as follows in Eq. 1:
Where N = total number of patients (80 for TCGA, 103 for ref. 82), k = number of co-mutations (2 for TCGA, 2 for ref. 82), x = total number of CYSLTR2 mutations (3 for TCGA, 2 for ref. 82) and y = total number of semaphorin/plexin mutations (4 for TCGA, 11 for ref. 82). We also independently confirmed this expression via Monte Carlo simulation. This value describes the probability of observing this enrichment of mutations by chance for the given number of mutations and patients. All other analysis of patient data presented in Supplementary Fig. 4 was performed using cBioPortal (RRID:SCR_014555)53,54.
Bulk RNA-sequencing analysis
HEK 293T cells were transfected as described above and harvested 36 h after transfection. RNA was isolated using miRNeasy (QIAGEN) and mRNA sequencing libraries were prepared according to manufacturer’s protocol using RNA using Illumina TruSeq Stranded mRNA library preparation kit (Illumina). Three independent biological replicates were obtained for transfections, lysates, and RNA preparation. Raw reads from bulk RNA-sequencing were trimmed with Trim Galore (RRID:SCR_011847) v0.4.4_dev (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) and quality-checked with FastQC v0.11.8 (RRID:SCR_014583) (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Trimmed reads were then aligned to the hg38 human reference genome with STAR aligner v2.5.3a118 (RRID:SCR_004463), and converted to gene counts with the analyzeRepeats.pl script in HOMER119 (RRID:SCR_010881). Differential expression: Gene counts were normalized and queried for differential expression using DESeq2 v1.30.0120 (RRID:SCR_015687). For each pairwise comparison, genes with fewer than 20 total raw counts across all samples were discarded prior to normalization, and genes with an absolute log2foldchange > 1 and an FDR-corrected p ≤ 0.05 were pulled as significant. Functional Enrichment: Genes were queried for treatment-specific functional enrichment using over-representation analysis (ORA) and gene set enrichment analysis (GSEA) in WebGestaltR v0.4.4121. Differentially expressed genes in each pairwise comparison were queried against the biological process ontology with ORA, while pairwise GSEAs were used to query pathway gene sets from KEGG (RRID:SCR_012773) (https://www.genome.jp/kegg/pathway.html), and MSigDB (http://www.gsea-msigdb.org/gsea/msigdb/collections.jsp).
Cancer cell drug adaptation assay
92.1 UM cells were treated with either 100 nM of the FAKi VS-4718 (provided by the Gutkind lab) or DMSO for a period of 21 days. Treatments and media were refreshed every 72 h. Cells were plated for each condition in triplicate, and three biological replicates of the complete experiment were performed. RNA was isolated using E.Z.N.A. Total RNA kit I (Omega) at intervals of 0, 3, 7, 14, and 21 days. mRNA sequencing libraries were prepared according to manufacturer’s protocol using RNA using Illumina TruSeq Stranded mRNA library preparation kit (Illumina).
Data availability
All information needed to reproduce the model are included in the Supplementary Information. RNAseq data from the experiments new to this manuscript are available from the Gene Expression Omnibus (GEO) with accession numbers GSE267152 and GSE267153. A Python notebook that includes all of the information needed to reproduce the model is available at https://github.com/StitesLab/GNAQ_model. A PDF version is also included to facilitate review, in case one does not have Python installed.
References
Perou, C. M. et al. Molecular portraits of human breast tumours. Nature 406, 747–752 (2000).
Shipp, M. A. et al. Diffuse large B-cell lymphoma outcome prediction by gene-expression profiling and supervised machine learning. Nat. Med. 8, 68–74 (2002).
Beroukhim, R. et al. Assessing the significance of chromosomal aberrations in cancer: methodology and application to glioma. Proc. Natl Acad. Sci. USA 104, 20007–20012 (2007).
Cancer Genome Atlas Research. The cancer genome atlas pan-cancer analysis project. Nat. Genet. 45, 1113–1120 (2013).
Alexandrov, L. B. et al. Signatures of mutational processes in human cancer. Nature 500, 415–421 (2013).
Marty, R. et al. MHC-I genotype restricts the oncogenic mutational landscape. Cell 171, 1272–1283.e1215 (2017).
Stites, E. C., Trampont, P. C., Ma, Z. & Ravichandran, K. S. Network analysis of oncogenic Ras activation in cancer. Science 318, 463–467 (2007).
Stites, E. C., Trampont, P. C., Haney, L. B., Walk, S. F. & Ravichandran, K. S. Cooperation between noncanonical Ras network mutations. Cell Rep. 10, 307–316 (2015).
Stites, E. C. & Shaw, A. S. Quantitative systems pharmacology analysis of KRAS G12C covalent inhibitors. CPT Pharmacomet. Syst. Pharm. 7, 342–351 (2018).
McFall, T. et al. A systems mechanism for KRAS mutant allele-specific responses to targeted therapy. Sci. Signal. 12, eaaw8288 (2019).
McFall, T., Schomburg, N. K., Rossman, K. L. & Stites, E. C. Discernment between candidate mechanisms for KRAS G13D colorectal cancer sensitivity to EGFR inhibitors. Cell Commun. Signal. 18, 179 (2020).
McFall, T. & Stites, E. C. Identification of RAS mutant biomarkers for EGFR inhibitor sensitivity using a systems biochemical approach. Cell Rep. 37, 110096 (2021).
McFall, T. et al. Co-targeting KRAS G12C and EGFR reduces both mutant and wild-type RAS-GTP. NPJ Precis. Oncol. 6, 86 (2022).
Wey, M., Lee, J., Jeong, S. S., Kim, J. & Heo, J. Kinetic mechanisms of mutation-dependent Harvey Ras activation and their relevance for the development of Costello syndrome. Biochemistry 52, 8465–8479 (2013).
Hunter, J. C. et al. Biochemical and structural analysis of common cancer-associated KRAS mutations. Mol. Cancer Res. 13, 1325–1335 (2015).
Donovan, S., Shannon, K. M. & Bollag, G. GTPase activating proteins: critical regulators of intracellular signaling. Biochim. Biophys. Acta 1602, 23–45 (2002).
Stites, E. C. & Ravichandran, K. S. Mechanistic modeling to investigate signaling by oncogenic Ras mutants. Wiley Interdiscip. Rev. Syst. Biol. Med. 4, 117–127 (2012).
Stites, E. C. & Ravichandran, K. S. Mathematical investigation of how oncogenic ras mutants promote ras signaling. Methods Mol. Biol. 880, 69–85 (2012).
Stites, E. C. Chemical kinetic mechanistic models to investigate cancer biology and impact cancer medicine. Phys. Biol. 10, 026004 (2013).
McFall, T. & Stites, E. C. A mechanism for the response of KRAS(G13D) expressing colorectal cancers to EGFR inhibitors. Mol. Cell Oncol. 7, 1701914 (2020).
Markevich, N. I. et al. Signal processing at the Ras circuit: what shapes Ras activation patterns? Syst. Biol. 1, 104–113 (2004).
Kiel, C. & Serrano, L. Cell type-specific importance of ras-c-raf complex association rate constants for MAPK signaling. Sci. Signal. 2, ra38 (2009).
Wolf, J., Dronov, S., Tobin, F. & Goryanin, I. The impact of the regulatory design on the response of epidermal growth factor receptor-mediated signal transduction towards oncogenic mutations. FEBS J. 274, 5505–5517 (2007).
Saucerman, J. J., Healy, S. N., Belik, M. E., Puglisi, J. L. & McCulloch, A. D. Proarrhythmic consequences of a KCNQ1 AKAP-binding domain mutation: computational models of whole cells and heterogeneous tissue. Circ. Res. 95, 1216–1224 (2004).
Lew, E. D., Furdui, C. M., Anderson, K. S. & Schlessinger, J. The precise sequence of FGF receptor autophosphorylation is kinetically driven and is disrupted by oncogenic mutations. Sci. Signal. 2, ra6 (2009).
Yeung, E. et al. Inference of multisite phosphorylation rate constants and their modulation by pathogenic mutations. Curr. Biol. 30, 877–882.e876 (2020).
Zewde, N. & Morikis, D. A computational model for the evaluation of complement system regulation under homeostasis, disease, and drug intervention. PLoS ONE 13, e0198644 (2018).
Kraikivski, P., Chen, K. C., Laomettachit, T., Murali, T. M. & Tyson, J. J. From START to FINISH: computational analysis of cell cycle control in budding yeast. NPJ Syst. Biol. Appl. 1, 15016 (2015).
Benary, U., Kofahl, B., Hecht, A. & Wolf, J. Modeling Wnt/beta-catenin target gene expression in APC and Wnt gradients under wild type and mutant conditions. Front Physiol. 4, 21 (2013).
Mitchell, S., Tsui, R., Tan, Z. C., Pack, A. & Hoffmann, A. The NF-kappaB multidimer system model: a knowledge base to explore diverse biological contexts. Sci. Signal. 16, eabo2838 (2023).
Basak, S. et al. A fourth IkappaB protein within the NF-kappaB signaling module. Cell 128, 369–381 (2007).
Chen, J., Yue, H. & Ouyang, Q. Correlation between oncogenic mutations and parameter sensitivity of the apoptosis pathway model. PLoS Comput. Biol. 10, e1003451 (2014).
Benedict, K. F. et al. Systems analysis of small signaling modules relevant to eight human diseases. Ann. Biomed. Eng. 39, 621–635 (2011).
Krantz, B. A., Dave, N., Komatsubara, K. M., Marr, B. P. & Carvajal, R. D. Uveal melanoma: epidemiology, etiology, and treatment of primary disease. Clin. Ophthalmol. 11, 279–289 (2017).
Takasaki, J. et al. A novel Galphaq/11-selective inhibitor. J. Biol. Chem. 279, 47438–47445 (2004).
Schrage, R. et al. The experimental power of FR900359 to study Gq-regulated biological processes. Nat. Commun. 6, 10156 (2015).
Onken, M. D. et al. Targeting nucleotide exchange to inhibit constitutively active G protein alpha subunits in cancer cells. Sci. Signal. 11, eaao6852 (2018).
Annala, S. et al. Direct targeting of Galpha(q) and Galpha(11) oncoproteins in cancer cells. Sci. Signal. 12, eaau5948 (2019).
Yang, J., Manson, D. K., Marr, B. P. & Carvajal, R. D. Treatment of uveal melanoma: where are we now? Ther. Adv. Med. Oncol. 10, 1758834018757175 (2018).
Jenkins, R. W. & Fisher, D. E. Treatment of advanced melanoma in 2020 and beyond. J. Investig. Dermatol. 141, 23–31 (2021).
Nathan, P. et al. Overall survival benefit with tebentafusp in metastatic uveal melanoma. N. Engl. J. Med. 385, 1196–1206 (2021).
Oldham, W. M. & Hamm, H. E. Heterotrimeric G protein activation by G-protein-coupled receptors. Nat. Rev. Mol. Cell Biol. 9, 60–71 (2008).
Robertson, A. G. et al. Integrative analysis identifies four molecular and clinical subsets in uveal melanoma. Cancer Cell 32, 204–220.e215 (2017).
Maziarz, M. et al. Atypical activation of the G protein Galpha(q) by the oncogenic mutation Q209P. J. Biol. Chem. 293, 19586–19599 (2018).
Linderman, J. J. Modeling of G-protein-coupled receptor signaling pathways. J. Biol. Chem. 284, 5427–5431 (2009).
Turcotte, M., Tang, W. & Ross, E. M. Coordinate regulation of G protein signaling via dynamic interactions of receptor and GAP. PLoS Comput. Biol. 4, e1000148 (2008).
Katanaev, V. L. & Chornomorets, M. Kinetic diversity in G-protein-coupled receptor signalling. Biochem. J. 401, 485–495 (2007).
Feng, X. et al. Hippo-independent activation of YAP by the GNAQ uveal melanoma oncogene through a trio-regulated rho GTPase signaling circuitry. Cancer Cell 25, 831–845 (2014).
Feng, X. et al. A platform of synthetic lethal gene interaction networks reveals that the GNAQ uveal melanoma oncogene controls the hippo pathway through FAK. Cancer Cell 35, 457–472.e455 (2019).
Carvajal, R. D. et al. Advances in the clinical management of uveal melanoma. Nat. Rev. Clin. Oncol. 20, 99–115 (2023).
Mukhopadhyay, S. & Ross, E. M. Rapid GTP binding and hydrolysis by G(q) promoted by receptor and GTPase-activating proteins. Proc. Natl Acad. Sci. USA 96, 9539–9544 (1999).
Wingler, L. M. et al. Angiotensin analogs with divergent bias stabilize distinct receptor conformations. Cell 176, 468–478.e411 (2019).
Gao, J. et al. Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci. Signal. 6, pl1 (2013).
Cerami, E. et al. The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov. 2, 401–404 (2012).
Patt, J. et al. An experimental strategy to probe Gq contribution to signal transduction in living cells. J. Biol. Chem. 296, 100472 (2021).
Ostrem, J. M., Peters, U., Sos, M. L., Wells, J. A. & Shokat, K. M. K-Ras(G12C) inhibitors allosterically control GTP affinity and effector interactions. Nature 503, 548–551 (2013).
Lito, P., Solomon, M., Li, L. S., Hansen, R. & Rosen, N. Allele-specific inhibitors inactivate mutant KRAS G12C by a trapping mechanism. Science 351, 604–608 (2016).
Skoulidis, F. et al. Sotorasib for lung cancers with KRAS p.G12C mutation. N. Engl. J. Med. 384, 2371–2381 (2021).
Ma, J., Weng, L., Bastian, B. C. & Chen, X. Functional characterization of uveal melanoma oncogenes. Oncogene 40, 806–820 (2021).
Chidiac, P. & Ross, E. M. Phospholipase C-beta1 directly accelerates GTP hydrolysis by Galphaq and acceleration is inhibited by Gbeta gamma subunits. J. Biol. Chem. 274, 19639–19643 (1999).
Moore, A. R. et al. Recurrent activating mutations of G-protein-coupled receptor CYSLTR2 in uveal melanoma. Nat. Genet. 48, 675–680 (2016).
Nell, R. J. et al. Involvement of mutant and wild-type CYSLTR2 in the development and progression of uveal nevi and melanoma. BMC Cancer 21, 164 (2021).
Johansson, P. et al. Deep sequencing of uveal melanoma identifies a recurrent mutation in PLCB4. Oncotarget 7, 4624–4631 (2016).
Moore, A. R. et al. GNA11 Q209L mouse model reveals RasGRP3 as an essential signaling node in uveal melanoma. Cell Rep. 22, 2455–2468 (2018).
Slater, K. et al. Evaluation of oncogenic cysteinyl leukotriene receptor 2 as a therapeutic target for uveal melanoma. Cancer Metastasis Rev. 37, 335–345 (2018).
Ceraudo, E. et al. Direct evidence that the GPCR CysLTR2 mutant causative of uveal melanoma is constitutively active with highly biased signaling. J. Biol. Chem. 296, 100163 (2021).
Shirley, M. D. et al. Sturge-Weber syndrome and port-wine stains caused by somatic mutation in GNAQ. N. Engl. J. Med. 368, 1971–1979 (2013).
Subramanian, A. et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl Acad. Sci. USA 102, 15545–15550 (2005).
Mootha, V. K. et al. PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat. Genet. 34, 267–273 (2003).
Cordenonsi, M. et al. The Hippo transducer TAZ confers cancer stem cell-related traits on breast cancer cells. Cell 147, 759–772 (2011).
Field, M. G. et al. PRAME as an Independent Biomarker for Metastasis in Uveal Melanoma. Clin. Cancer Res. 22, 1234–1242 (2016).
Gene Ontology, C. The Gene Ontology resource: enriching a GOld mine. Nucleic Acids Res. 49, D325–D334 (2021).
Ashburner, M. et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat. Genet. 25, 25–29 (2000).
Mi, H., Muruganujan, A., Ebert, D., Huang, X. & Thomas, P. D. PANTHER version 14: more genomes, a new PANTHER GO-slim and improvements in enrichment analysis tools. Nucleic Acids Res. 47, D419–D426 (2019).
Jassal, B. et al. The reactome pathway knowledgebase. Nucleic Acids Res. 48, D498–D503 (2020).
Worzfeld, T. & Offermanns, S. Semaphorins and plexins as therapeutic targets. Nat. Rev. Drug Discov. 13, 603–621 (2014).
Sakurai, A. et al. Semaphorin 3E initiates antiangiogenic signaling through plexin D1 by regulating Arf6 and R-Ras. Mol. Cell Biol. 30, 3086–3098 (2010).
Ng, T. et al. Class 3 semaphorin mediates dendrite growth in adult newborn neurons through Cdk5/FAK pathway. PLoS One 8, e65572 (2013).
Bechara, A. et al. FAK-MAPK-dependent adhesion disassembly downstream of L1 contributes to semaphorin3A-induced collapse. EMBO J. 27, 1549–1562 (2008).
Field, M. G. et al. BAP1 loss is associated with DNA methylomic repatterning in highly aggressive class 2 uveal melanomas. Clin. Cancer Res. 25, 5663–5673 (2019).
Urtatiz, O., Haage, A., Tanentzapf, G. & Van Raamsdonk, C. D. Crosstalk with keratinocytes causes GNAQ oncogene specificity in melanoma. Elife 10, e71825 (2021).
Johansson, P. A. et al. Whole genome landscapes of uveal melanoma show an ultraviolet radiation signature in iris tumours. Nat. Commun. 11, 2408 (2020).
Adzhubei, I. A. et al. A method and server for predicting damaging missense mutations. Nat. Methods 7, 248–249 (2010).
Argast, G. M. et al. Plexin B1 is repressed by oncogenic B-Raf signaling and functions as a tumor suppressor in melanoma cells. Oncogene 28, 2697–2709 (2009).
Chakraborty, G., Kumar, S., Mishra, R., Patil, T. V. & Kundu, G. C. Semaphorin 3A suppresses tumor growth and metastasis in mice melanoma model. PLoS ONE 7, e33633 (2012).
Stevens, L. et al. Plexin B1 suppresses c-Met in melanoma: a role for plexin B1 as a tumor-suppressor protein through regulation of c-Met. J. Investig. Dermatol. 130, 1636–1645 (2010).
Paradis, J. S. et al. Synthetic lethal screens reveal cotargeting FAK and MEK as a multimodal precision therapy for GNAQ-driven uveal melanoma. Clin. Cancer Res. 27, 3190–3200 (2021).
Wong, O. G. et al. Plexin-B1 mutations in prostate cancer. Proc. Natl Acad. Sci. USA 104, 19040–19045 (2007).
Balakrishnan, A. et al. Molecular profiling of the “plexinome” in melanoma and pancreatic cancer. Hum. Mutat. 30, 1167–1174 (2009).
Martin, M. et al. Exome sequencing identifies recurrent somatic mutations in EIF1AX and SF3B1 in uveal melanoma with disomy 3. Nat. Genet. 45, 933–936 (2013).
Diederichs, S. et al. The dark matter of the cancer genome: aberrations in regulatory elements, untranslated regions, splice sites, non-coding RNA and synonymous mutations. EMBO Mol. Med. 8, 442–457 (2016).
Jung, H., Lee, K. S. & Choi, J. K. Comprehensive characterisation of intronic mis-splicing mutations in human cancers. Oncogene 40, 1347–1361 (2021).
Lyon, A. M. & Tesmer, J. J. Structural insights into phospholipase C-beta function. Mol. Pharm. 84, 488–500 (2013).
Navaratnarajah, P., Gershenson, A. & Ross, E. M. The binding of activated Galpha(q) to phospholipase C-beta exhibits anomalous affinity. J. Biol. Chem. 292, 16787–16801 (2017).
Gutenkunst, R. N. et al. Universally sloppy parameter sensitivities in systems biology models. PLoS Comput. Biol. 3, 1871–1878 (2007).
Bender, R. J. & Mac Gabhann, F. Dysregulation of the vascular endothelial growth factor and semaphorin ligand-receptor families in prostate cancer metastasis. BMC Syst. Biol. 9, 55 (2015).
Zhang, L., Polyansky, A. & Buck, M. Modeling transmembrane domain dimers/trimers of plexin receptors: implications for mechanisms of signal transmission across the membrane. PLoS ONE 10, e0121513 (2015).
Calmelet, C., Madamanchi, A. & Zutter, M. Multiscale coupled modeling of Plexin-D1 and notch signaling in retinal sprouting angiogenesis. J. Coupled Syst. Multiscale Dyn. 5, 1–17 (2017).
Tanaka, T. et al. Hybrid in vitro/in silico analysis of low-affinity protein-protein interactions that regulate signal transduction by Sema6D. Protein Sci. 31, e4452 (2022).
Sakumura, Y., Tsukada, Y., Yamamoto, N. & Ishii, S. A molecular model for axon guidance based on cross talk between rho GTPases. Biophys. J. 89, 812–822 (2005).
Hu, Q. & Shokat, K. M. Disease-causing mutations in the G protein Galphas subvert the roles of GDP and GTP. Cell 173, 1254–1264.e1211 (2018).
Chen, Q. et al. Structures of rhodopsin in complex with G-protein-coupled receptor kinase 1. Nature 595, 600–605 (2021).
Masuho, I. et al. A global map of G protein signaling regulation by RGS proteins. Cell 183, 503–521.e519 (2020).
Inoue, A. et al. Illuminating G-protein-coupling selectivity of GPCRs. Cell 177, 1933–1947.e1925 (2019).
Gremer, L. et al. Germline KRAS mutations cause aberrant biochemical and physical properties leading to developmental disorders. Hum. Mutat. 32, 33–43 (2011).
Kiel, C. & Serrano, L. Structure-energy-based predictions and network modelling of RASopathy and cancer missense mutations. Mol. Syst. Biol. 10, 727 (2014).
Gillies, T. E. et al. Oncogenic mutant RAS signaling activity is rescaled by the ERK/MAPK pathway. Mol. Syst. Biol. 16, e9518 (2020).
Moller, I. et al. Activating cysteinyl leukotriene receptor 2 (CYSLTR2) mutations in blue nevi. Mod. Pathol. 30, 350–356 (2017).
Goto, K., Pissaloux, D., Paindavoine, S., Tirode, F. & de la Fouchardiere, A. CYSLTR2-mutant Cutaneous melanocytic neoplasms frequently simulate “pigmented epithelioid melanocytoma,” expanding the morphologic spectrum of blue tumors: a clinicopathologic study of 7 cases. Am. J. Surg. Pathol. 43, 1368–1376 (2019).
van de Nes, J. et al. Targeted next generation sequencing reveals unique mutation profile of primary melanocytic tumors of the central nervous system. J. Neurooncol. 127, 435–444 (2016).
Kusters-Vandevelde, H. V. N. et al. Whole-exome sequencing of a meningeal melanocytic tumour reveals activating CYSLTR2 and EIF1AX hotspot mutations and similarities to uveal melanoma. Brain Tumor Pathol. 35, 127–130 (2018).
van de Nes, J. A. P. et al. Activating CYSLTR2 and PLCB4 mutations in primary leptomeningeal melanocytic tumors. J. Investig. Dermatol. 137, 2033–2035 (2017).
Zhou, J. et al. Somatic mutations of GNA11 and GNAQ in CTNNB1-mutant aldosterone-producing adenomas presenting in puberty, pregnancy or menopause. Nat. Genet. 53, 1360–1372 (2021).
Akin-Bali, D. F. Bioinformatics analysis of GNAQ, GNA11, BAP1, SF3B1,SRSF2, EIF1AX, PLCB4, and CYSLTR2 genes and their role in the pathogenesis of uveal melanoma. Ophthalmic Genet. 42, 732–743 (2021).
Yang, J. H. et al. A white-box machine learning approach for revealing antibiotic mechanisms of action. Cell 177, 1649–1661.e1649 (2019).
AlQuraishi, M. & Sorger, P. K. Differentiable biology: using deep learning for biophysics-based and data-driven modeling of molecular mechanisms. Nat. Methods 18, 1169–1180 (2021).
Herman, J. & Usher, W. SALib: an open-source Python library for sensitivity analysis. J. Open Source Softw. 2, 97 (2017).
Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).
Heinz, S. et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol. Cell 38, 576–589 (2010).
Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).
Liao, Y., Wang, J., Jaehnig, E. J., Shi, Z. & Zhang, B. WebGestalt 2019: gene set analysis toolkit with revamped UIs and APIs. Nucleic Acids Res. 47, W199–W205 (2019).
Acknowledgements
We would like to thank Nasun Hah and the members of the Next Generation Sequencing (NGS) Core at Salk, Ling Huang, and the Razavi Newman Integrative Genomics and Bioinformatics (IGC) Core, and Amy Cao. We would also like to thank all members of the Stites lab for helpful comments and advice, in particular Thomas McFall and Jianfeng Huang for offering helpful experimental design suggestions, and Gaurav Mendiratta for helpful discussions on statistical analysis. This work was supported by NIH T32CA009370 (M.T.), NIH DP2AT011327 (E.C.S.), the Melanoma Research Alliance Young Investigator Award (E.C.S.), the Conrad Prebys Foundation (E.C.S.), by the Razavi Newman Integrative Genomics and Bioinformatics and the NGS Core Facility of the Salk Institute with funding from NIH-NCI CCSG: P30 014195, the Chapman Foundation, and the Helmsley Charitable Trust.
Author information
Authors and Affiliations
Contributions
Conceptualization: M.Tr., E.C.S.; Performed computational analysis: M.Tr., K.L.; M.To., M.B. Performed experimental analysis: M.Tr., N.K., K.A.; Writing—original draft: M.Tr., E.C.S.; Writing—review and editing: M.Tr., K.A., N.A., M.To., M.B., J.S.G., E.C.S.; Resources: M.B., J.S.G., E.C.S.; Supervision: J.S.G., E.C.S.; Funding acquisition: E.C.S.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Trogdon, M., Abbott, K., Arang, N. et al. Systems modeling of oncogenic G-protein and GPCR signaling reveals unexpected differences in downstream pathway activation. npj Syst Biol Appl 10, 75 (2024). https://doi.org/10.1038/s41540-024-00400-1
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41540-024-00400-1
