Neuroblastoma signalling models unveil combination therapies targeting feedback-mediated resistance

Abstract

Very high risk neuroblastoma is characterised by increased MAPK signalling, and targeting MAPK signalling is a promising therapeutic strategy. We used a deeply characterised panel of neuroblastoma cell lines and found that the sensitivity to MEK inhibitors varied drastically between these cell lines. By generating quantitative perturbation data and mathematical modelling, we determined potential resistance mechanisms. We found that negative feedbacks within MAPK signalling and via the IGF receptor mediate re-activation of MAPK signalling upon treatment in resistant cell lines. By using cell-line specific models, we predict that combinations of MEK inhibitors with RAF or IGFR inhibitors can overcome resistance, and tested these predictions experimentally. In addition, phospho-proteomic profiling confirmed the cell-specific feedback effects and synergy of MEK and IGFR targeted treatment. Our study shows that a quantitative understanding of signalling and feedback mechanisms facilitated by models can help to develop and optimise therapeutic strategies. Our findings should be considered for the planning of future clinical trials introducing MEKi in the treatment of neuroblastoma.

Fig 1. Mutations are insufficient to explain sensitivity variations to RAS/PI3K drugs in neuroblastoma cell line panel.

A: Oncoprint of 9 neuroblastoma cell lines for RAS/p53/PI3K related genes along with MYCN and TERT mRNA expression. Bold font indicates MYCN-amplified cell lines. B: Relative IC50 of the same 9 neuroblastoma cell lines as in A for drugs targeting the PI3K and MAPK pathways (n = 2). C: Viability concentration curves for the MEK inhibitor AZD6244 on the neuroblastoma cell line panel along with the calculated IC50 (intersection with dotted line). Points represent measurements (n = 2).

Fig 2. Neuroblastoma cell lines show heterogeneous responses to signalling perturbations.

A: Outline of the perturbation experiments. A panel of cell lines was treated with growth factors and small molecule inhibitors, and the resulting effect on selected phosphoproteins was measured using multiplexed bead-based ELISAs. B: Graphical representation of the perturbation scheme on a literature signalling network. Blue and red contour highlights ligand stimulation and kinase inhibition, respectively; yellow filling shows measured phosphoproteins. C: Perturbation data obtained from applying all combinations of 4 ligands or BSA control and 7 inhibitors or DMSO control to 6 neuroblastoma cell lines. Each measurement is normalised by the BSA+DMSO control of the corresponding cell line and represents at least 2 biological replicates. Readouts are phospho-proteins p-MEK1^S217/S221, p-p38^T180/Y182, p-ERK1^T202/Y204, p-cJUN^S63, p-AKT^S473 and p-S6K^T389. D: Global non-mechanistic analysis of the perturbation data presented in C: top first two components of a principal component analysis and bottom hierarchical clustering. Colour scale corresponds to the IC50 for AZD6244 treatment (see also Fig 1C).

Fig 3. Receptor expression and topology variations explain the heterogeneity in perturbation response.

A: Starting from a literature-derived network, a model was fitted for each cell line (Initial model fit) and extended following suggestions from the model (Model extensions and refit). Those models with different network structures were then harmonised by fixing the inhibition parameters to a consensus value (Fixed inhibitor parameters) to make the parameters directly comparable (Parameter comparison). B: Model residuals before and after model extension and harmonisation. The black line represents the number of data points, which is equal to the expected mean of the error if the model explains all the data. C: Cell-line-specific network extensions (dashed arrows) relative to the literature network. Colour of the extended link was matched to cell line colour if required in only one cell line model and black otherwise. D: Model paths from the receptors to the first measured downstream node and correlation with the corresponding receptor expression. The colours correspond to the value of the path scaled by the maximum absolute value of that path between all cell lines. E: Model paths between non-receptor perturbed nodes and measured nodes for routes present in at least 2 cell lines. Colour scale is the same as in D. Cells are ordered from left to right from most sensitive to most resistant to the MEK inhibitor AZD6244. Due to the absence of ASK1 basal activity in IMR32, ASK1->p38 and ASK1->MEK represent in this cell line NGF->ASK1->p38 and NGF->ASK1->MEK, respectively.

Fig 4. AZD6244 resistant cell lines have strong feedback control of MAPK signalling.

A: Mean pMEK log2-fold change relative to control after AZD6244 treatment in 7 neuroblastoma cell lines measured with bead-based ELISAs. Error bars represent 95% confidence interval. B: Measurement of 6 phosphoproteins (columns) after perturbation of N206, IMR32 and KELLY by either EGF (KELLY, N206) or NGF (IMR32) (together referred to as GF), IGF1, or control BSA in combination with Sorafenib (RAFi), AZD6244 (MEKi) or control DMSO. Values are expressed in log2-fold change to BSA+DMSO control. C: Starting model and S6K→IGF1 receptor extension for the high pMEK responder cell lines. D: (top panel) Model residuals for N206, IMR32 and KELLY models with (black) or without (blue) an S6K→IGF1 receptor feedback link and corresponding p-value(χ² test with df = 1). (bottom panel) Parameter values of the high pMEK responder models including the S6K→IGF1 receptor link.

Fig 5. Phosphoproteomics analysis reveals important variations in the response to combination treatment.

A-B: Venn diagrams showing the overlap in differentially regulated phosphosites A: between IMR32 and N206 or B: between treatments for each cell line. C: Phosphopeptides synergistically altered by MEKi+IGFRi combination (black outline) when compared to the sum of individual inhibitor treatments. AKT, mTOR or P70S6K bona fide targets (bold font) and putative targets (italic font; top 5 predicted kinases by PhosphoNET Kinase Predictor www.phosphonet.ca) are indicated. D: Kinase substrate enrichment score using PhosphoSitePlus annotations. E: Log-fold change to DMSO for RAF/MAPK and MYCN phosphopeptides. C-E: Black outline highlights significant changes in activity (limma moderated t-test, FDR<5%) F-H: Relative levels compared to control of the total proteins levels for MYCN (F) and CCND1 (H) measured with mass spectrometry and MYCN measured with Western blot (G).

Fig 6. AZD6244 resistant cell lines can be sensitised with combined inhibition with the IGFR inhibitor AEW541 or the RAF inhibitor LY3009120.

A: Model-inferred targeting strategy of dual inhibition. B: Growth inhibition measurements for various combinations of the MEK inhibitor AZD6244 with the RAF inhibitor LY3009120 or the IGFR inhibitor AEW541. Values over 100 indicate cell death. n = 2. C: Bliss synergy corresponding to the measurements in B. D: left: Model predictions of pERK and pMEK activity for MEK inhibition alone and in combination with inhibition of upstream kinase RAF or downstream kinase ERK for KELLY and IMR32. Values are log-fold changes to IGF1 condition with inhibitor strength set to -1. D: right: pERK and pMEK plex measurements in KELLY and IMR32 after 90min treatment of the MEK inhibitor AZD6244 in combination with either DMSO, SCH772984 (ERKi, 10μM) or LY3009120 (RAFi, 5μM) in cells grown with 10% FCS. Values are log-fold change to FCS medium condition.