Abstract
Val600Glu (V600E) mutation is the most common BRAF mutation detected in thyroid cancer. Hence, recent research efforts have been performed trying to explore several inhibitors of the V600E mutation-containing BRAF kinase as potential therapeutic options in thyroid cancer refractory to standard interventions. Among them, vemurafenib is a selective BRAF inhibitor approved by Food and Drug Administration for clinical practice. Unfortunately, vemurafenib often displays limited efficacy in poorly differentiated and anaplastic thyroid carcinomas probably because of intrinsic and/or acquired resistance mechanisms. In this view, cancer stem cells (CSCs) may represent a possible mechanism of resistance to vemurafenib, due to their self-renewal and chemo resistance properties.
We present a computational framework to suggest new potential targets to overcome drug resistance. It has been validated with an in vitro model based upon a spheroid-forming method able to isolate thyroid CSCs that may mimic resistance to vemurafenib. Indeed, vemurafenib did not inhibit cell proliferation of BRAF V600E thyroid CSCs, but rather stimulated cell proliferation along with a paradoxical over-activation of ERK and AKT pathways. The computational model identified a fundamental role of mitogen-activated protein kinase 8 (MAP3K8), a serine/threonine kinase expressed in thyroid CSCs, in mediating this drug resistance. To confirm model prediction, we set a suitable in vitro experiment revealing that the treatment with MAP3K8 inhibitor restored the effect of vemurafenib in terms of both DNA fragmentation and poly (ADP-ribose) polymerase cleavage (apoptosis) in thyroid CSCs. Moreover, MAP3K8 expression levels may be a useful marker to predict the response to vemurafenib.
The model is available in GitHub repository visiting the following URL: https://github.com/francescopappalardo/MAP3K8-Thyroid-Spheres-V-3.0.
Supplementary data are available at Bioinformatics online.
1 Introduction
Activating mutation in the BRAF serine/threonine kinase represents the most common genetic alteration in thyroid cancer, occurring in approximately 45% of papillary thyroid cancer (PTC) and in a lower proportion of poorly differentiated thyroid cancer (PDTC) and anaplastic thyroid cancer (ATC) (Begum et al., 2004; Nikiforova et al., 2003; Xing, 2005).
Over 95% of all BRAF mutations detected in thyroid cancer is a thymine to adenine transversion at exon 15 nucleotide 1799 (T1799A) of the BRAF gene leading to substitution of valine by glutamic acid at residue 600 of the protein chain (V600E) (Trovisco et al., 2006). This alteration results in a constitutively active BRAF molecule with sustained kinase activity that, in turn, promotes chronic stimulation of mitogen-activated protein kinase (MAPK) pathway, thereby resulting in increased phosphorylation of downstream targets, including ERK kinase (Nikiforov and Nikiforova, 2011). This altered signaling results in increased cell proliferation, survival and tumor progression, in a growth factor independent manner (Knauf and Fagin, 2009). Also, other BRAF mutations in thyroid cancer are rarely reported (Carta et al., 2006; Ciampi et al., 2005; Hou et al., 2007).
Given the crucial role in thyroid cancer, BRAF V600E represents a potential therapeutic target for treatment of patients with advanced thyroid cancer refractory to standard approaches. In fact, although surgery, radioiodine and thyroid-stimulating hormone-suppressive therapy result in improved outcomes for most patients with well-differentiated thyroid carcinomas (papillary and follicular), effective options for the treatment of advanced disease including radioiodine-resistant and metastatic differentiated thyroid cancer, poorly differentiated and anaplastic thyroid carcinomas remain limited.
An emerging therapeutic strategy for these cancers is targeted therapy with small molecule kinase inhibitors (KIs). Several KIs have been explored for the treatment of advanced thyroid tumors, but just some of them have been approved for the use in clinical practice. To date, the US Food and Drug Administration (FDA) approved sorafenib and lenvatinib for the treatment of advanced differentiated thyroid carcinoma. In contrast, KIs did not display any significant effect in patients with poorly differentiated and ATC. Among them, also a BRAF V600E inhibitor vemurafenib (PLX4032, approved by FDA for treatment of metastatic melanoma) did not provide any satisfactory result as a single agent in patients with poorly differentiated and anaplastic thyroid carcinomas. In a small cohort of seven patients with ATC, only two patients achieved durable responses to vemurafenib therapy while four patients had progression of disease (Hyman et al., 2015). This result suggests that poorly differentiated and ATC cells display mechanisms tumor relapse due to intrinsic/acquired resistance to vemurafenib.
A possible explanation of relapse and resistance to conventional therapies in poorly differentiated and anaplastic thyroid disease has been provided by the so-called cancer stem cells (CSCs). This is a relatively small cell sub-population within tumor mass with stem cell like properties and the ability to grow as non-adherent spheroids (thyrospheres) and to sustain self-renewal in vitro ( Ahn et al., 2014; Gianì et al., 2015; Malaguarnera et al., 2011).
Thyrospheres isolated from thyroid cancer tissue and continuous cell lines are characterized by their ability to sustain tumor growth, promote metastasis and confer resistance to chemotherapeutic agents ( Giuffrida et al., 2016; Li et al., 2013; Liotti et al., 2017; Reeb et al., 2014; Todaro et al., 2010).
To elucidate mechanisms of resistance to vemurafenib and define strategies to more efficiently target BRAF V600E, we established a reliable spheroid-forming assay to identify and enrich thyroid CSCs, then tested the effect of vemurafenib on these cells.
Computational modeling and in silico simulations represent now-a-days essential resources for the analysis and the investigation of controversial and unknown mechanisms of resistance to conventional drug therapy and are able to predict possible new interventions to overcome this issue still unresolved. Existing works about several computational approaches have been aimed to improve knowledge about resistance in drug targets (Pan et al., 2016; Sun and Hu, 2018). We developed a computational model based on a system of differential equations (in particular, we used ordinary differential equations, (ODEs)) to analyze the entire protein network of MAPK and PI3K/AKT pathways, in which BRAF V600E mutation is strongly involved. In silico simulations, coupled with graph search algorithms, suggested a prominent role of tumor progression locus 2 (MAP3K8) kinase, also known as COT, in resistance of 8505C-CSCs to vemurafenib. The model prediction about MAP3K8 was then confirmed by in vitro experiments. Moreover, MAP3K8 might represent a potential novel prognostic marker and therapeutic target for advanced thyroid cancer treatment to further investigate.
2 Materials and methods
2.1 Cell cultures, thyrospheres formation and self-renewal assay
ATC cells HTH74 (BRAF wild-type) and SW1736 (BRAF V600E positive) were provided by Dr. N. E. Heldin (Uppsala, Sweden). The 8505-C (BRAF V600E positive) cell line was purchased from European Animal Cell Culture (Salisbury, UK). The follicular cancer cell line WRO, expressing wild-type BRAF gene, was provided by A. Fusco (Naples, Italy). BRAF V600E mutation in thyroid cancer cell lines was confirmed by PCR amplification and sequencing, as previously reported (Frasca et al., 2008). These cell lines were cultured as monolayers at 37°C with 5% CO2 under humidified conditions in RPMI-1640 medium (Sigma) supplemented with 10% (v/v) heat inactivated fetal bovine serum (FBS; Thermo), 2 mM L-glutamine and 1% penicillin-streptomycin (Sigma).
To generate thyrospheres, cells were trypsinized from adherent cultures, washed three times in large volumes of calcium-magnesium-free PBS, then seeded as single cells at density 1–5 x 104 cells/ml in serum-free DMEM/F12 medium supplemented with 20 ng/ml epidermal growth factor ((EGF); Sigma) and 10 ng/ml basic fibroblast growth factor ((bFGF); Thermo) and B-27 supplement (1: 50; Thermo) on low-attachment flasks (Falcon). Half of medium volume was replaced every 3 days. After 7–10 days, spherical-like floating groups of cells greater than 50 μm were scored and considered as spheres.
For self-renewal assay, spheres were collected by gentle centrifugation and dissociated enzymatically with 0.05% trypsin/EDTA. Single cells were plated at density 1 x 104 cells/ml on six-well low-attachment plates to generate secondary thyrospheres. The efficiency of sphere formation was expressed as percentage after 7–10 days of observation.
2.2 Computational model
The model is built on a system of ODEs able to simulate both the expected and the emergent behavior within the protein network.
We developed the in silico model including pathways and their associated protein networks, with the molecular entities and interactions involved in CSCs resistance to vemurafenib. In this view, we tailored resistance to vemurafenib by focusing on the mechanisms underlying the enhanced ERK/AKT phosphorylation in response to the drug.
2.2.1 Identification of entities and biochemical interactions
The starting point to select the entities and the biochemical reactions needed to simulate the ERK dynamics in 8505 C CSCs was to find the appropriate signaling pathways cascade. To this aim, we instantiated the model using as main source the Kyoto Encyclopedia of Genes and Genomes PATHWAY Database (KEGG) (Kanehisa et al., 2004). The following steps were carried out on KEGG database: (i) obtain all the entities that participate in the biochemical reactions involved in tumorigenesis of thyroid cancer considering the protein kinase cascades that sequentially contain EGF-Sos-Ras-Raf1-Mek1-ERK and EGFR-PI3K-PIP3-AKT; (ii) retrieve all the nodes and the biochemical interactions from the KEGG map owning cross-talks with EGF-Sos-Ras-Raf1-Mek1-ERK and EGFR-PI3K-PIP3-AKT cascades; (iii) include into the computational model the interactions between the two specific pathways, ‘Kegg reference map ko04010’ (for MAPK signaling pathway), ‘Kegg reference map hsa04151’ (for PI3K/AKT). The protein kinase cascades involved in EGFR activation, PIK3CA function, AKT modulation and RAF1 signaling were also included. Moreover, RAS signaling pathway (Kegg reference: ko04014) and mTOR signaling pathway (Kegg reference: ko041150) were included within the model to simulate mechanisms in terms of possible escapement of apoptosis and alternative processes that may induce Vemurafenib resistance phenomena through the activation of additional signaling pathways.
Serine/threonine-protein kinase BRAF represents the main step in the transduction of mitogenic signals from the cell membrane to the nucleus. BRAF protein is strongly linked both to MAPK and PI3K/AKT pathways that are mainly involved in tumor progression. Hence, we focused on BRAF protein and its role in regulating the signaling pathways involved in cell proliferation, differentiation and apoptosis. Moreover, we performed a comparative modeling on how wild type BRAF gene and its V600E mutant variant may influence the entire biochemical system.
The model tailored for 8505 C CSCs was built with 75 entities, including protein kinases [i.e. RAC serine/threonine-protein kinase (AKT)], growth factors [i.e. EGF], membrane receptors in their different status of activation (i.e. pEGFR in its activated form), small molecules, including protein KIs at different concentrations (i.e. vemurafenib) and 85 biochemical reactions.
To simulate the dynamics of the system, we used the COmplex PAthway Simulator (COPASI) pathway software tool (Mendes et al., 2009), commonly used in systems biology. All the ODEs and the relative kinetic functions are reported in Supplementary Figure S1. The source code of the model is publicly available via GitHub server, in the link specified in the availability section. Moreover, we provided an ‘how to use’ guide to facilitate the usage of the model (Supplementary Document S2).
Besides the classic Henri–Michaelis Menten law, which includes the majority of enzyme-catalyzed biochemical reactions in biology (Bajzer and Strehler, 2012), we used a slightly modified version of this law that is more suitable for our modeling purposes. It is able to analyze the ratio between the rate of enzymatic reactions and their affinity for the substrates, taking in account the efficiency of the modifier for each catalysis that participates in the kinetic law as a further substrate.
As shown in the Equation (2), this rule develops towards an ‘inverted’ direction: if the modifier (substrate1) is absent or at low concentration, the reaction could take place. Otherwise, if the modifier is in a normal or constitutive activated status, the biochemical reaction will not take place.
In order to find ever long-distance crosstalk in the graph defined by the signaling network under study, we applied depth first search approach (DFS) (Cormen et al., 2002). DFS is an algorithm that is commonly employed as a search algorithm for data structures based on trees and graphs. The name of the algorithm inherits its meaning from the type of search: it proceeds to visit nodes that are far from the root before starting the visit of closer nodes. This property is a key feature of DFS. In comparison to other graph search algorithms such as breadth-first search (BFS) that starts visiting nodes from the root, DFS starts visiting nodes from the leaves. In the event that the goal coincides with the identification of nodes more likely closer to root, BFS would represent a better choice. However, if the target node is close to a leaf, as in this case, we would use DFS. This choice allows to maximize the possibility to detect long-distance nodes connected to the root.
Specifically, when applied to graphs, DFS identifies trees that are sub-graphs. In these sub-graphs, it is possible to distinguish forward edges that are directed from one node of the graph to one of its descendants and back edges that, instead, are directed from one node to one of its ancestors (inverted effect).
We applied DFS to the MAPK pathway, selecting ERK as the starting point i.e. the root, to automatically discover nodes that participate in the modulation of ERK activity.
Among all the returned sub-graphs (i.e. sub-pathways) that DFS produced in output, we considered the one depicted in Figure 1 where it is possible to identify two main arms. The first one originates from RAS signaling and the other one comes from TNF signaling. Specifically, the MAP3K8 node influences ERK through MEK signaling. We modeled this specific sub-pathway in order to analyze the effect of MAP3K8 in the ERK modulation.
One of the sub-pathway returned by DFS (MAPK, ERK) call. DFS algorithm received in input the MAPK signaling pathway and the ERK node as root. Among the sub-graphs returned by the search, the one depicted in the figure shows an interesting node (MAP3K8) that modulates ERK activity through MEK arm bypassing RAS-BRAF path
2.2.2 Model parameters setting
To reproduce the experimental conditions of cell environment, we assumed that all the molecular entities examined had a homogenous concentration overall the cell volume, ranging from ≈ 10 nM/l to ≈10 μM/l (Legewie et al., 2008). For mutated species, in particular BRAF, we set the concentration to 100 nM/l, because this protein is constitutively activated and overexpressed during the entire time course of the simulation; for growth factors entities that occupy the culture medium [i.e. Fibroblast growth factor (FGF)] the concentration was set to 1 x 104 nM/l; for all the remaining molecules to 10 nM/l.
The number of enzymatic reactions catalyzed per second and the Henri–Michaelis–Menten constant that measures the affinity between the enzyme and the substrate were set to the same value i.e. 0.1. This is because experimental data do not indicate any predominant velocity in the considered biochemical reactions.
To reproduce drug resistance phenomenon, we introduced the Mutated-BRAF species at the concentration of 100 nM/l that is the same concentration of the initial Inactive-BRAF. We then removed all the biochemical pathways not influencing Mutated-BRAF species (i.e. the activation of BRAF by Rap1) and replaced Active- BRAF with Mutated-BRAF species in the switch of Mek activation. Moreover, new elements were subsequently inserted to reproduce the inhibitory effect of vemurafenib. To this purpose, we introduced vemurafenib species at a concentration of 1000 nmol/l that corresponds to the same one used in vitro in 8505C-CSCs. To better define the reaction framework of vemurafenib, we also included its degradation rate and its inhibitory effect against BRAF Mutated-species. Indeed, vemurafenib half-life in vivo is about 57 h (Zhang et al., 2017) and this parameter was used to rate the associated mass action law and simulate its decay.
2.3 Reagents, real time PCR, cell growth, western blot and ELISA quantification
All details regarding used reagents, real time PCR analysis, Western Blot preparation and ELISA quantification of DNA fragmentation are available in Supplementary Document S1.
3 Results
3.1 Enrichment and identification of a CSC sub-population from different thyroid cancer cells lines
We first identified and characterized CSCs from thyroid cancer cell lines based on their ability to form spheres during stem cells culture conditions. Then their self-renewal potential and the expression of CSC-associated genes were evaluated.
To this purpose single cell suspension from four PDTC cell lines –8505 C, SW1736, WRO and HTH74 ––were cultured at low density into low-attachment plates in serum-free stem cell conditioned culture medium. After 6–10 days, floating spheroids with different morphological phenotypes were observed in each culture (Fig. 2 A). All our cell lines contained sphere-forming cells, which in turn gave rise to secondary spheres after dissociation. Moreover, the spheres had been serially passaged up to six generations, providing the evidence of their self-renewal capability in vitro. To assess the phenotype of cells within the spheres, we measured the expression of CSCs markers OCT3/4, ABCG2, HEY-1 and HES-1 by real-time PCR. When compared to monolayer parental cells, thyrospheres expressed significantly higher levels of mRNA for all of the genes examined (Fig. 2B).
These data indicate that a sub-population of cancer stem-like cells with self-renewal ability is present in thyroid cancer cell lines, and it is enriched in sphere culture conditions.
Isolation and characterization of CSCs from different PDTC cell lines. (A) Representative phase-contrast microscopy images of floating spheroid isolated from four PDTC cell lines. Cells were placed in serum-free culture medium containing EGF and bFGF in non-adherent culture flasks to form spheres as indicated in methods. (B) qRT-PCR analysis of stem cell markers. Data were normalized for GAPDH expression and related to parental cells. Data shown represent mean + /– SEM from at least three independent experiments. Statistical significance was determined by t-test (* = P < 0.05; ** = P < 0.01)
3.2 Different effect of vemurafenib on MAPK signaling and cell growth in thyroid cancer cell monolayer and sphere counterparts
To evaluate the effect of BRAF inhibitor vemurafenib, we first compared the effect of 1 μM vemurafenib on ERK activation, in both wild-type and BRAF V600E mutated monolayer cultures and thyrospheres, respectively. As expected (Fig. 3 A), vemurafenib inhibited ERK phosphorylation in cell monolayer carrying BRAF V600E mutation (8505 C and SW1736), whereas it failed to inhibit ERK signaling in wild-type BRAF cell monolayer (WRO and HTH74), confirming the selective effect of vemurafenib for BRAF V600E. Unexpectedly, in mutant BRAF V600E spheres (derived from 8505 C and SW1736 cells) vemurafenib only transiently inhibited ERK phosphorylation. Moreover, vemurafenib had no effect on ERK phosphorylation in BRAF wild type cells, both monolayer and spheres cells (WRO and HTH74).
Effect of vemurafenib on BRAF mutant and wild type CSC and monolayer cells counterparts. (A) Effects of vemurafenib on ERK1/2 phosphorylation in parental (left panel) and spheres (right panel) thyroid cancer cells. The indicated cell lines were incubated with 1 μM for 2 h. Lysates were then separated by SDS-PAGE and immunoblotted with the indicated antibodies. Representative blots and densitometric analysis of three independent experiments are shown. (B) Effect of vemurafenib on 8505C (expressing mutant BRAF) cell number in both spheres (dotted bars) and parental monolayer cells (gray bars). Representative blots and densitometric analysis of three independent experiments are shown (*P < 0.05 w.r.t. control; **P < 0.01 w.r.t. control; ***P < 0.001 w.r.t. control). Data shown represent mean + /– SEM from at least three independent experiments
Thus, we compared the effect of vemurafenib on growth inhibition in mutant BRAF V600E both monolayer 8505 C cells and spheres. Consistent with ERK phosphorylation inhibition by vemurafenib observed by Western blot in adherent 8505 C cells, treatment with vemurafenib caused a significant reduction on cells number (by 32% see Fig. 3B). In contrast, vemurafenib did not inhibit cell proliferation in 8505C-spheres, but an increased cell proliferation was rather observed compared to control (Fig. 3B).
These data indicate that, compared to BRAF V600E mutated cell monolayer, CSC-enriched spheres respond to vemurafenib displaying a paradoxical cell proliferation.
3.3 Effect of vemurafenib on cell signaling in 8505C-CSCs
In order to identify the key molecular mediators involved in the 8505C-CSC resistance to vemurafenib, we next investigated the kinetics of main steps of signaling pathways.
To this end, time course experiments (4 h) exploring the levels of both ERK and AKT activation, were performed in 8505C-CSCs treated with vemurafenib. Western blot analysis showed that, at variance with 8505 C cell monolayer, the inhibition of ERK was transient, with a paradoxical re-activation beginning 30 min after addition of the drug in 8505C-CSCs. Interestingly, western blot analysis of AKT phosphorylation revealed a concomitant progressive increase in AKT activation (Fig. 4 ). Unlike 8505C-CSCs, vemurafenib significantly inhibited ERK phosphorylation in standard monolayer culture of 8505 C. No significant effect was observed in phosphorylation of AKT.
Time course experiment of the effect of vemurafenib on signaling in BRAF mutated cells. 8505C-monolayer cells (left panel) and 8505C-spheres (right panel) were treated with 1 μM vemurafenib and collected at the indicated times. Lysates were then separated by SDS-PAGE and immunoblotted with the indicated antibodies. Representative blots and densitometric analysis of three independent experiments are shown (*P < 0.05 w.r.t. control; **P < 0.01 w.r.t. control; ***P < 0.001 w.r.t. control)
3.4 Model predictions on the mechanisms of drug resistance in 8505C-CSCs
The developed model allows to investigate alternative signaling pathways, even over long distances, applying exploring methodologies based on DFS in MAPK and PI3K/AKT signaling pathways.
We identified other pathways that do not directly influence ERK dynamics. In particular, we found that TNF signaling pathway that is divided in two branches, stimulates both ERK and AKT protein kinases; ERK pathway is activated by TNFR1, while AKT by TNFR2. This bifurcation contains fundamental nodes, not considered before, to stimulate ERK and AKT phosphorylation. ERK activation is allowed through MAP3K8 contribution. AKT activation is triggered through the intervention of the TRAF complex and then, at a second stage, through PI3K activation. At this stage, the DFS carried on our in silico model of 8505C-CSCs allowed to predict that the concomitant inhibition of both MAP3K8 and BRAF V600E is able to consistently reduce ERK phosphorylation. This prediction was confirmed by in vitro experiments with selective inhibitors (see the next sub-section). In silico predictions about the inhibitor combination protocol are shown in Figure 5 where the computational framework was applied to simulate and predict three different cases. The first (panel A) shows the behavior of pERK in 8505C-CSCs thyrospheres under BRAF mutation condition. The model correctly reproduced what observed in vitro i.e. an uncontrolled increasing of pERK that leads to a cellular mechanism of resistance to vemurafenib. The contribution of the species in activating pERK i.e. bRafMutated, simulates the constitutively activation of BRAF V600E. Moreover, Mek and Raf1 species enhance ERK phosphorylation.
In silico prediction of pERK dynamics on 8505C-CSCs after MAP3K8 + vemurafenib administration. Panel A shows pERK concentration at time 4 h while panel B (linear scale) and panel C (logarithmic scale) depict its dynamics over time. ERK behavior under BRAF knock-out condition when vemurafenib is administered at 1 µM concentration is simulated setting the PLX4032 species value to 1000 nM in the in silico model. These results show an initial p-ERK de-activation after BRAF arm knock-out and then a significant increase of the resistance mechanism against vemurafenib and the ineffectiveness of the treatment with conventional KIs. When MAP3K8 inhibitor is administered at 10 µM concentration along with vemurafenib (at 1 µM concentration) the in silico model predicts a significant decrease of pERK (more than half of reduction, in comparison with the vemurafenib only administration). Inhibition of MAP3K8 was simulated setting MAP3K8: NF-kB_Inactive species value to 0 in the in silico model
Panel B describes pERK activation when 1 µM vemurafenib is added. An initial pERK inactivation after BRAF arm inhibition and then the appearance of the ERK over-activation in response to vemurafenib renders ineffective the BRAF inhibitor action. PLX4032 (vemurafenib) species affects bRafMutated species. This simulates the inhibitory activity of the treatment over BRAF V600E. The bRafMutated neighborhood (Mek and Raf1 species) continues to provide its contribution to pERK activation, even if vemurafenib is administered. Last panel (panel C) shows the model predictions when MAP3K8 10 µM inhibitor is added along with 1 µM vemurafenib. In this case, pERK activation is significantly decreased, suggesting an anti-resistance mechanism elicited by the inhibition of the target predicted by the model i.e. MAP3K8. MAP3K8 species activates Mek1 leading to ERK phosphorylation. MAP3K8 knock-out strongly contributes to pERK de-activation, predicting resistance avoidance.
The entire modeled signaling pathway is reported in Supplementary Figure S2.
We performed dynamical analysis to assess the significance of trends through stochastic simulations. The stochastic model represents the number of particles of each chemical species and use a reaction probability density function (PDF) to describe the timing of reaction events. Gillespie (Gillespie, 2007), developed a Monte Carlo simulation algorithm, known as the stochastic simulation algorithm (SSA) first reaction method, that simulates the stochastic dynamics of the system by sampling this PDF. As a single simulation using SSA reports just one grasp of a probabilistic representation, providing limited information, we ran 50 simulations for the untreated case (vemurafenib only) and 50 simulations for the treated one (vemurafenib + MAP3K8 inhibitor). This allowed to analyze both dynamic evolution and the significance of trends. For the untreated case, we obtained a mean pERK concentration of 16.76 nM/ml (SD + /–1.5); for the treated case, we obtained a mean pERK concentration of 2.924 nM/ml (SD + /–2.438).
Moreover, a two tailed Wilcoxon signed rank test for matched pairs was executed to assess the efficacy of the treatment on pERK dynamics. The null hypothesis H0 is that measurements are drawn from the same population, so possible differences come purely from chance. The test confirmed the rejection of the null hypothesis (P < 0.0001), confirming a strong statistical difference between the two groups (treated versus untreated). Figure 6 depicts stochastic simulations results of pERK dynamics for both untreated and treated cases.
In silico stochastic simulation of pERK dynamics on 8505C-CSCs in both untreated (vemurafenib only administration) and treated (MAP3K8 + vemurafenib administration) cases. Fifty different simulations were performed for each case. For the untreated case, the mean value of pERK concentration was of 16.76 nM/ml (SD + /–1.5). For the treated case, the mean value of pERK concentration was of 2.924 nM/ml (SD + /–2.438)
To analyze in deeper details the contribution of pathways species to pERK dynamics we further performed sensitivity analysis. Generally speaking, sensitivity analysis provides a mean to measure how much a selected model variable (the effect) changes when a selected parameter (the cause) is changed. In particular, the sensitivity coefficients describe how fast species concentrations at steady state change when the values of the kinetic parameters are slightly increased.
Through COPASI simulation platform, we set a sensitivity analysis directed to investigate the effects of variations in species concentrations to pERK dynamics, in steady state condition, for both the untreated and treated cases. Figure 7 depicts the results of the sensitivity analysis. Panel A and B refer to the treated case, while Panel C and D refer to untreated case. Panel A shows the sensitivity coefficients related to species that, when perturbed, lead to an increase of pERK steady state concentration. In particular, the pathway arms related to the contribution of IKK complex, TRADD complex, TAB complex, NIK, AKT, MAP3K8, PI3K and Mek species play an important role in the modulation of BRAF mutation when constitutively activated in 8505C-CSCs. Panel B reports the sensitivity coefficients related to species that, when perturbed, lead to a decrease of pERK steady state concentration. Specifically, the sensitivity analysis reveals that MAP3K8 inhibition is the main cause of pERK reduction. Panel C depicts the species that, in the untreated case, contribute to pERK increase. It is worth mentioning that in this case small perturbation of MAP3K8 leads to pERK increase, confirming its prominent role. Finally, panel D shows how perturbations of RAS, RAF, PTEN, EGF and FGF arms still make a non-negligible mean contribution to pERK decrease. In other words, these pathway arms are essential to ERK activation machinery. This is in good agreement with previous literature results, supporting the model robustness.
Sensitivity analysis results. In the treated case, the sensitivity analysis shows how particular species lead to an increase (panel A) or to a decrease (panel B) of pERK. In particular IKK complex, TRADD complex, TAB complex, NIK, AKT, MAP3K8, PI3K and Mek modulate an increase of pERK while MAP3K8 inhibition is the main cause of pERK reduction. In the untreated case, panel C reveal sensitivity coefficients of species that affect pERK increase (MAP3K8 along with NIK, TAB complex, IKK complex and TRAD complex) while panel D depicts the scenario in which pathway arms containing RAS, RAF, PTEN, EGF and FGF are essential to ERK activation machinery
3.5 MAP3K8 is up-regulated in CSCs and inhibition of its kinase activity increases sensitivity to vemurafenib in thyroid CSCs harboring BRAF mutation
As a first step toward understanding the role of MAP3K8 kinase in resistance of 8505C-CSCs to vemurafenib, we compared the expression levels of MAP3K8 between the 8505C-CSCs and monolayer 8505 C cancer cells counterpart by qRT-PCR. Interestingly, 8505C-CSCs showed higher levels of MAP3K8 expression than 8505 C monolayer parental cells (Fig. 8 A).
MAP3K8 kinase expression and function in 8505C-Spheres. (A) qRT-PCR of MAP3K8 expression in 8505C-spheres (dotted bar) and 8505C-monolayer (gray bar) cells. Data, normalized for GAPDH expression and compared to parental cells, represent mean + /– SEM from at least three independent experiments. Statistical significance was determined by t-test (** = P < 0.01). (B) Effect of MAP3K8 inhibitor on mutant BRAF 8505C-spheres proliferation in response to vemurafenib. Spheres formation was measured 6 days after addition of the indicated concentration of MAP3K8 inhibitor in the absence (gray bars) or the presence of 1 μM vemurafenib (dotted bars). Data are expressed as percent change (mean +/– SEM) compared to untreated cells. (C, D) Effect of MAP3K8 and BRAF V600E inhibitor leads on 8505C-sphereapoptosis. Western blot analysis of cleaved PARP (c) and colorimetric assay of apoptotic DNA fragmentation (d) were evaluated in 8505C-spheres treated with either vemurafenib (gray bar), MAP3K8 inhibitor (hatched bar) or both (dotted bar). Cells were harvested at 24 h after drug treatment and evaluated as indicated in methods. (E) 8505C-spheres were treated with 1 μM vemurafenib alone or in combination with 10 μM MAP3K8 inhibitor, and collected at the indicated times. Lysates were then subjected to SDS-PAGE and immunoblotted with the indicated antibodies. Representative blots and densitometric analysis of three independent experiments are shown
Next, we evaluated the effect of a selective inhibitor of MAP3K8 on 8505C-CSCs response to vemurafenib in terms of cell growth and apoptosis by measuring both the cleavage of poly (ADP-ribose) polymerase (PARP) and DNA fragmentation analysis.
To this end, we exposed 8505C-CSCs to increasing doses of the MAP3K8 inhibitor with or without 1 µM vemurafenib for 7 days. MAP3K8 inhibitor significantly restored in a dose-dependent manner the effect of vemurafenib on sphere formation in 8505C-CSCs harboring the BRAF V600E mutation (Fig. 8B).
According to growth inhibition, western blot analysis revealed a substantial PARP cleavage in 8505C-CSCs upon combined treatment, indicative of apoptosis induction (Fig. 8C). More specifically, treatment of MAP3K8 inhibitor used as single agent did induce PARP cleavage; however, it was less marked as compared to combined treatment. As expected, vemurafenib alone failed to induce cleavage of PARP. Similar results were observed in DNA fragmentation assay (Fig. 8D).
Finally, combined treatment of vemurafenib and MAP3K8 inhibitor resulted in effective suppression of ERK rebound and AKT over-activation compared to vemurafenib alone, providing a rationale for the observed synergy in both growth and apoptotic assays (Fig. 8E). We also confirmed our results using an additional thyroid cancer cell lines harboring BRAF mutation, SW1736 ( Supplementary Fig. S3). Moreover, we explored whether the presence of the MAP3K8 KI increases sensitivity to vemurafenib of both monolayer 8505 C and SW1736 cancer cells, harboring BRAF V600E mutation. To this end, cells were incubated with increasing doses of MAP3K8 inhibitor in presence or absence of 1 μM of the selective BRAF inhibitor vemurafenib. As shown in Supplementary Figure S4, MAP3K8 inhibitor significantly enhanced the inhibitory effect of vemurafenib on cell growth in thyroid cancer cell lines harboring BRAF V600E. However, we observed that the MAP3K8 KI alone displayed a more pronounced effect in spheres than monolayer cancer cell counterpart. Indeed, the IC50 values were 7.8 μM in 8505C-spheres and 6.8 μM in SW1736-spheres (see Fig. 8B and Supplementary Fig. S3B, respectively), while the IC50 value >10 μM on both monolayer cancer cell lines.
4 Discussion and conclusions
The present results indicate that a small sub-population of cancer cells, CSCs, may be involved in the failure of therapies targeting BRAF V600E in advanced thyroid cancer. These data have been obtained in CSCs derived from a panel of well-established thyroid cancer cell lines by a sphere-formation assay. This cell system provided a reliable model for further experiments on KI resistance. The stem cell identity of these cancer spheres is supported by their extensive self-renewal in vitro and overexpression of several stemness-related genes, such as OCT-4, ABCG-2, HES-1 and HEY-1, compared to monolayer parental cells.
Here, we find a paradoxical stimulatory effect of BRAF-inhibitor vemurafenib on mutant BRAF V600E spheres (8505C-spheres) at variance with their parental monolayer cultures (8505C-monolayers). Indeed, same drug dose, while inhibiting cell growth and ERK activation in mutant BRAF V600E monolayer cultures, it stimulated sphere proliferation with a markedly different effect on kinetics of both ERK and AKT activation. In 8505C-spheres vemurafenib induced only a transient inhibition of ERK phosphorylation with a rapid and paradoxical re-activation starting at 30 min. In the same experimental context, a progressive increase in AKT activation was also observed. In accordance with the theory of the role of CSCs in chemo-resistance, these data suggest a resistance phenomenon to vemurafenib in a sub-population of these thyroid cancer cells.
Importantly, the in silico approach allowed us to identify MAP3K8 protein, a serine/threonine kinase expressed in 8505C-CSCs, as possible mechanism of resistance. Here, we developed a computational model with the aim to investigate alternative signaling transduction pathways over long distances that can affect possible crosstalk links between MAPK and PI3K/AKT.
Applying DFS strategy, we tracked alternative pathways able to explain the escapement from ERK activation. Under BRAF inhibition, the DFS identified long distance entities that, in turn, may affect the dynamics of ERK activation. By a computational view, this method identifies forward edges, which point from a node of the graph to one of its descendants (for example TNF for ERK and TRAF complex for AKT) and back edges, which point from a node to one of its ancestors (inverted effect). By this search, we found that MAP3K8 (MAP3K8) contributes to ERK activation, suggesting that this kinase may be targeted to overcome drug resistance.
In vitro experiments indicated that in mutant BRAF V600E 8505C-CSCs the inhibitory activity of vemurafenib is restored by the addiction of the selective MAP3K8 inhibitor. Hence, these data were in agreement with and confirmed our prediction obtained by simulation modeling.
MAP3K8 kinase was shown to be an important mediator in inflammation during both innate and adaptive immune response ( Arthur and Ley, 2013; Lee et al., 2015c; Miyoshi et al., 1991). Activation of MAP3K8 involves downstream signaling pathways such as MEK/ERK, JNK, p38 MAPK and NF-κB not only involved in inflammatory response but also in the production of pro-inflammatory cytokines such as TNFα, IL-1β and IFNγ (Gantke et al., 2011; Martel and Rousseau, 2014; Pattison et al., 2016; Vougioukalaki et al., 2011). Recently, accumulating evidences suggest that MAP3K8 may play a role in cell transformation, tumor growth and metastatic progression in several human cancers ( Jeong et al., 2011; Kim et al., 2014; Lee et al., 2015b; Lee et al., 2016; Sourvinos et al., 1999). In this view, a recent report indicates that that higher expression of MAP3K8 correlates with both the BRAF V600E mutation and tumor recurrence in PTC (Lee et al., 2015a).
Our data indicate for the first time that MAP3K8 is up-regulated in 8505C-CSCs compared to 8505 C cell monolayer. Moreover, treatment with MAP3K8 inhibitor alone induced both marked DNA fragmentation and PARP cleavage, indicative of apoptosis induction in 8505C-CSCs. Furthermore, the findings obtained by our computational model are in agreement with those described in ovarian cancer cells, where MAP3K8 levels are positively related to the ERK phosphorylation levels (Gruosso et al., 2015). Moreover, other previous data, obtained in melanoma cells, further confirm our findings. Indeed, in this cell model MAP3K8 expression is associated with de novo resistance in BRAF V600E cultured cell lines and acquired resistance to BRAF inhibitors (Johannessen et al., 2010),
Taken together, MAP3K8 appears not only drive vemurafenib resistance in CSCs but might represent a potential novel prognostic marker and therapeutic targets for advanced thyroid cancer treatment. Indeed, MAP3K8 expression levels may be a clinically useful biomarker to predict the response to vemurafenib monotherapy in BRAF mutant tumors as well as may also help guide the selection of MAP3K8 combination therapy in these cancers.
Moreover, these data are suggestive that MAP3K8 may support the development of stem cell capabilities in the thyroid cancer cells and, also be associated with tumor progression and metastasis. Thus, future study will be directed to investigate the role of MAP3K8 in thyroid cancer.
In conclusion, our data highlight a distinct mechanism by which mutant BRAF V600E thyroid cancer escape from vemurafenib inhibitory activity. Our current data reveal that MAP3K8 contributes to vemurafenib resistance of mutant BRAF V600E thyroid CSCs.
The further role of MAP3K8 as biomarker in predicting response in thyroid cancer under vemurafenib treatment could be better assessed through in vivo experiments that will be performed in due course.
Moreover, with a limited effort, the in silico model can be adapted to other diseases that share the same signaling pathway.
Conflict of Interest: none declared.
References
Author notes
The authors wish it to be known that, in their opinion, Fiorenza Gianì and Giulia Russo authors should be regarded as Joint First Authors.
