Neuroblastoma signalling models unveil combination therapies targeting feedback-mediated resistance

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.


Introduction
Neuroblastoma is the most common and devastating extracranial childhood solid tumour, accounting for    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 Figure 1C  When we applied hierarchical clustering on the cell line panel, SKNSH was clustered separately, suggesting 125 that it has a very atypical response to the perturbations, with a generally very high response to all ligands, 126 and an especially strong response to PDGF ( Figure 2D bottom). This atypical status of SKNSH is also 127 5 . CC-BY-NC 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted June 14, 2021. present in the mRNA expression, with a PCA on the most variables genes or on the genes in the GO term "signal transduction" separating it from the other cell lines. Interestingly, CHP212 also separated from the other cell line in a PCA based on gene expression data, but not when considering the response to 130 the perturbations. When grouping cells by MEK inhibitor sensitivity, we noticed that simple multivariate 131 analysis by PCA does not separate cells into groups that correspond to sensitive or resistant cells ( Figure 2D

134
Signalling models highlight differential feedback regulation of MEK

135
To get further, more mechanistic, insights into potential resistance mechanisms, we used the perturbation 136 data to parameterise signalling models. We applied our previously developed method that has been derived 137 from Modular Response Analysis (MRA, implemented as R package STASNet, Dorel et al (2018)  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-linespecific 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. 6 . CC-BY-NC 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted June 14, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 signalling network models to each cell line. This modelling procedure requires a literature network and the perturbation data as input, and then estimates response coefficients corresponding to link strengths using a 140 maximum likelihood estimate (see Figure 3A, first step). By using the statistical framework of the likelihood 141 ratio test, the modelling procedure then allows to test if any extension of the literature network is required 142 to describe the data (see Figure 3A, second step). To compare parameters between cell lines, it is essential 143 to harmonise parameters between all cells that can practically not be identified alone, i.e. parameters for 144 inhibitors (see Figure 3A, third step). This finally yields a parameter map that allows to compare signalling 145 strength between cell lines (see Figure 3A, final step).

146
When starting with a canonical literature network (see Materials and Methods), we obtained reason-147 able fits for 4 of the 6 cell lines, as judged by the sum of weighted squared residuals that is in the 148 order of number of data points ( Figure 3B, red bars), and the normal distribution of residuals (Sup-149 plementary Figure 11). When we systematically tested if extensions of the network improve the fit us-150 ing a likelihood ratio test, we found that significant improvements were still possible for most cell lines. 151 We therefore performed successive rounds of extensions for each cell line independently ( Figure 3A    which might indicate that NGF signalling is mediated mostly via NTRK1 in those cell lines. The parameters 175 for IGF-induced signals correlated with IGF1R or IGF2R for MEK and AKT, respectively, indicating that 176 both receptors mediate IGF1 signalling independently. Interestingly, the parameters for the pathway from 177 EGF to MEK did not correlate with EGFR expression, but they do for EGF to AKT, which might suggest 178 that differences in adaptor protein expression shape routing into downstream signalling in the various cell CC-BY-NC 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted June 14, 2021. ; https://doi.org/10. 1101/2021 In contrast to the receptor-associated parameters, the strength of intra-cellular kinase paths are less variable, and most paths are comparable between cell lines ( Figure 3E). The most prominent exception 194 is the negative feedback in MAPK signalling from ERK to RAF. When compared to the other cell lines, 195 this feedback appears to be 3 to 4 times stronger in KELLY and IMR32, which are two cell lines that are 196 highly resistant to AZD6244. A strong RAF-mediated feedback is a known resistance mechanism against  Figure 4A).

216
To more precisely dissect the feedback wiring, we generated additional focused perturbation data for those 217 cells with high feedback (KELLY, IMR32 and N206) to MEK inhibition. We stimulated cells with different 218 growth factors (IGF and NGF or EGF), and blocked MAPK signalling with MEK and RAF inhibitors, and 219 subsequently monitored six phosphoproteins ( Figure 4B). Subsequently, we used this data to parameterise a 220 focused MRA model that additionally either contained or did not contain the only receptor-mediated feedback 221 found in the first modelling round from S6K→IGF1 ( Figure 3C and Figure 4A). Inclusion of the IGF receptor-222 mediated feedback led to a significantly better fit of the data for N206 and KELLY (χ 2 p<0.05), but did 223 not improve the IMR32 model ( Figure 4C and D). Interestingly, the S6K→IGF1→RAF→MEK feedback is 224 stronger in the N206 models, but the pathway-intrinsic feedback (ERK→RAF→MEK) is stronger in KELLY 225 ( Figure 4D). This highlights that all these cells display negative feedback regulation, but the strengths of 226 the two layers of feedbacks are different between cell lines. (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  . CC-BY-NC 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made When we investigated the phosphorylation of components of the MAPK pathway more closely, we found 262 many RAF negative feedback/crosstalk sites to be down-regulated after MEK inhibition (BRAF: T401, S750, 263 T753; RAF1: S29, S642, S259) in both cell lines ( Figure 5E). MEK1 S222/S226 phosphorylation is increased 264 and pERK S204 decreased in both cell lines after MEK inhibition, in line with corresponding measurements 265 using bead-based ELISAs. Among those down-regulated phosphosites that were only significant in the 266 combination in N206 we detected many MYCN-phosphosites, notably MYCN S62, which is regulated by 267 MAPK via CDK1 (?). Interestingly, this loss of S62 phosphorylated MYCN is associated with reduced 268 MYCN levels ( Figure 5F). This downregulation was observed in IMR32 and N206 cells upon single inhibition 269 (IGFRi for N206 and MEKi for both cell lines), but only in N206 cells an even stronger downregulation could 270 be observed upon double inhibition ( Figure 5F). We confirmed these effects in Western blots for IMR32 and 271 N206 cells (Figure 5G), and also found downregulation of MYCN upon IGFRi as well as MEKi treatment 272 but no synergetic decrease after the combination treatment ( Figure 5G). Another interesting protein that is 273 regulated synergistically in N206 is Cyclin D1 (Figure 5H), a protein that is involved in cell cycle progression 274 and whose loss likely mediates MYCN loss. It should be noted that only 5 proteins (PHGDH, DERL1, 275 AMPD3, ARHGEF16 and CCND1) were found differentially affected with an FDR < 10%, highlighting that 276 on this time scale phospho-protein changes dominated.

277
Taken together, the proteomics data is coherent with the model that MAPK signalling in N206 is con- in N206 but little in KELLY or IMR32 ( Figure 6B).

291
When trying to overcome the model-derived strong ERK-RAF feedback found in all three cell lines with 292 a combination of MEK and RAF inhibition we only found a synergistic effect for two of the three cell lines 293 (N206 and KELLY), whereas IMR32 remained resistant and no synergy could be detected. We hypothesised 294 that this observed resistance in IMR32 might be either because the vertical inhibition by MEKi and RAFi was 295 molecularly not effective or that IMR32 might no longer depend on ERK signalling for survival and growth.

296
To distinguish the former from the latter we decided to compare model simulation and measurements for 297 perturbation effects of selected inhibitor combinations on pMEK and pERK in IMR32 and KELLY cells.

298
Based on the model simulations, in both cell lines the vertical inhibition of MEK + RAF inhibitor was 299 predicted to suppress MAPK signalling much stronger than MEK inhibitor alone or in combination with an 300 ERK inhibitor. Moreover, the suppressive effect was predicted to be even more profound in IMR32 than in 301 KELLY ( Figure 6C top). We then measured the effect on pMEK and pERK of MEK inhibitor alone and 302 in combination with the RAF inhibitor LY3009120 or ERK inhibitor SCH772984 (Figure 6C bottom). The 303 measurements qualitatively supported the model simulations showing that RAF inhibitor suppressed MEK 304 feedback activation by AZD6244, and that this suppression is stronger in IMR32. Addition of the ERK 305 inhibitor neither suppressed this feedback activation nor could it decrease ERK phosphorylation more than 306 RAF inhibition, as also predicted by the model. This suggests that in agreement with the model simulations 307 the combination of RAFi and MEKi is most effective in IMR32 to effectively suppress ERK activation and 308 feedback-mediated re-activation. However, since the growth is least affected by this combination IMR32 309 seems not to depend on ERK activity.

310
In the end, we identified 2 combinations effective at low drug concentrations against the MEK-inhibitor 311 10 . CC-BY-NC 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  Figure 5: Phosphoproteomics analysis reveals important variations in the response to combination treatment Venn diagram showing the overlap in differentially regulated phosphosites A between IMR32 and N206 or B between treatments for each cell line. C Phoshpopeptides synergistically altered by MEK+IGFR combination. Black outline highlights where the change in the combination is significantly different to the sum of the individual changes (limma moderated t-test, FDR<5%). D Kinase substrate enrichment score using phosphositeplus annotations ; black outline highlights significant changes in activity for a given condition (limma moderated t-test, FDR<5%) E Log-fold change to DMSO for RAF/MAPK and MYCN phosphopeptides ; black outline shows significantly altered phosphosites per condition (limma moderated t-test, FDR<5%). F-H Relative levels compared to control of the total proteins levels, MYCN measured with mass spectrometry (F), Western blot (G), and CCND1 measured with mass spectrometry (H). 11 . CC-BY-NC 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted June 14, 2021. ; https://doi.org/10.1101/2021.06.14.448322 doi: bioRxiv preprint resistant cell lines KELLY and N206. As both KELLY and N206 have strong multi-layered feedbacks (Fig-312 ure 4D), we reasoned that a combination of IGFRi, RAFi and MEKi might be even more efficient as it 313 targets both feedbacks, irrespective of their individual strength. We thus tested the effect of a combination 314 of AEW541, AZD6244 and LY3009120 and observed a >80% reduction in viability of both KELLY and N206 315 already at moderate concentration of all three drugs (300nM of AEW541, 50nM of LY3009120 and 500nM of 316 AZD6244) making it a potential therapeutic option ( Figure 6D and Supp_data_fig1_drug_sensitivity_fig5_synergies.zip).

317
12 . CC-BY-NC 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made   Figure 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 assessment by model simulations on pERK activity of 3 AZD6244 resistant neuroblastoma cell lines under various levels of MEK inhibition and IGFR or RAF inhibition B. Corresponding growth inhibition measurements using the specified inhibitors. n=2. C. top: 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. C. bottom: 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. D. Viability of the cell lines for selected concentrations of dual and triple inhibitor treatments targeting MEK, RAF and IGFR.

13
. CC-BY-NC 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted June 14, 2021. ; https://doi.org/10.1101/2021.06.14.448322 doi: bioRxiv preprint

Discussion
Neuroblastoma is a complex disease with distinct subtypes that display radically different outcomes, ranging 319 from spontaneous regression in low-risk groups to only 50% survival of patients in the high risk neuroblastoma 320 group. Mutations in RAS/MAPK signalling are a hallmark of high risk neuroblastoma, and also define a 321 subgroup of patients with ultra-high-risk neuroblastoma and an even worse survival. Therefore targeted 322 treatment might be a valid strategy to treat those patients. However, response to MEK inhibitors are very 323 variable, and it is thus important to understand mechanisms of resistance and how to circumvent these.

324
In this work, we explored how a more quantitative understanding of signalling can be used to design 325 combinatorial treatments to counteract drug resistance. We used a panel of deeply profiled cell lines rep-326 resenting high risk neuroblastoma and showed that the response to MEK inhibitors is variable, with some 327 cell lines responding at low doses in the nM range, whereas others are highly resistant. By using signalling

335
Our work highlights that systematic perturbation data are a powerful source to probe intracellular sig-336 nalling pathways. The connectivity of signalling pathways implies that minor quantitative alterations of the 337 network can lead to many changes in response, not all of which alter the phenotype. In this work, we saw 338 that multivariate analysis of the perturbation data alone was not fruitful to separate cell lines with respect 339 to their drug sensitivity. In contrast, integration of data by models highlighted that variations of only a few 340 links is enough to explain the differences between those cell lines. Modelling was therefore key to integrate 341 the data and to unveil feedback loops as potential sources of resistance.

342
In our work we used a maximum likelihood version of MRA, but there are multiple other methods 343 that might be suited to reconstruct semi-quantitative signalling networks from perturbation data. Oates  In this work, we showed that some neuroblastoma cell lines possess two major layers of feedback in 365 MAPK signalling. One of these feedbacks is pathway-intrinsic (from ERK to RAF) and one is a feedback to 366 the IGF receptor. Interestingly, different cell lines show different relative strength of feedbacks from ERK 367 to RAF and IGFR, and simulations show that those require different strategies for vertical inhibition. For 368 the cell line KELLY, modelling unveiled an extremely strong negative feedback from ERK to RAF. This 369 suggests that a combination of MEK and RAF inhibitor will be more potent than a combination of MEK 370 14 . CC-BY-NC 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted June 14, 2021. ; https://doi.org/10.1101/2021.06.14.448322 doi: bioRxiv preprint and IGFR inhibitor. In contrast, in the cell line N206, both feedbacks have similar strength, suggesting that both combinations might be potent. In line with these predictions, experiments showed that in KELLY 372 indeed the combination of MEK and RAF inhibitors is much more potent to reduce growth compared to the 373 combination of MEK and IGFR. In contrast, in N206 both combinations reduce growth.

374
Our phospho-proteomics analysis shows that the combination of MEK and IFGR also has different effects 375 in the two cell lines: Whereas it shows clearly synergistic effects of the combination in N206, there is no 376 sign of synergy in IMR32. By aggregating the phosphoproteome to kinase activities using kinase enrichment 377 scores, one can also get insight into the re-wireing of signalling after perturbation. In our case, it clearly 378 shows how the re-activation of RAF after MEK inhibition is inhibited by the treatment with IGFR inhibitors.

463
The p-RSK1 (S380) readout was discarded because of a low dynamic range.

464
The same procedure and analytes were used for the other perturbation assays in this paper. Refer to the 465 main text for the exact inhibitors and concentrations used for each experiment.

467
The model for each cell line was fitted separately from the corresponding perturbation data with the cre-ateModel function from the R package STASNet ((Dorel et al , 2018), https://github.com/molsysbio/ STASNet/releases/tag/Dorel2020). STASNet implements the variation of Modular Response Analysis (MRA) described in Klinger et al (2013) and Dorel et al (2018) that implements a dual effect of inhibitors as both a negative stimulus and a disruption of signal propagation. Under the hypothesis of pseudo-steady-state and locally linear dependencies between nodes, MRA models the response to a perturbation as where R ij is the global response of node j after perturbation of node i, r k ij is the local response of node j 468 after perturbation of node i taking into account the effect of inhibition of node k, and S ik is the sensitivity 469 of node i to perturbation k. For the proteomics and phosphoproteomics cells were grown to confluency for 3 days in full medium and 494 treated with AEW541 10µM and/or AZD6244 10µM or control DMSO for 4h.

495
17 . CC-BY-NC 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted June 14, 2021. ; https://doi.org/10.1101/2021.06.14.448322 doi: bioRxiv preprint a combination of LysC (Wako) and Trypsin (Promega) using the the single-pot, solid-phase-enhanced sample 497 preparation (?). For each sample, an equal amount of peptide was then chemically labelled with TMTpro 498 reagents (?). Samples were randomly assigned to one of the first 15 TMT channels, while the 16th channel was 499 composed of a superset of all the samples to allow multi-plex normalisation. Equal amounts of the labelling 500 reactions were combined in two TMT16 plexes, desalted via SepPak columns (Waters) and fractionated via 501 high-pH fractionation (?) on a 96 minutes gradient from 3 to 55% acetonitrile in 5 mM ammonium formate, 502 each fraction collected for 1 minute then combined into 24 fractions. From each fraction, an aliquot was 503 used to measure the total proteome while the remaining peptides were combined into 12 fractions and used 504 as input for an immobilised metal affinity chromatography using an Agilent Bravo system. For the total 505 proteome analysis, peptides were on-line fractionated on a multi-step gradient from 0 to 55% acetonitrile in 506 0.1% formic acid prior injection in a QExactive HF-x mass spectrometer. Samples were acquired using a data 507 dependent acquisition strategy with MS1 scans from 350 to 1500 m/z at a resolution of 60 000 (measured 508 at 200 m/z), maximum injection time (IT) of 10 ms and an automatic gain control (AGC) target value of 509 3 × 10 6 . The top 20 most intense precursor ions with charges from +2 to +6 were selected for fragmentation 510 with an isolation window of 0.7 m/z. Fragmentation was done in an HCD cell with a normalised collision 511 energy of 30% and analysed in the detector with a resolution of 45 000 (200 m/z), AGC target value of 10 5 , 512 maximum IT of 86 ms. We used the same parameters for phosphoproteome analysis with the exception of 513 MS2 maximum IT that was set to 240 ms.

514
The acquired raw files were analysed using MaxQuant v1.6.10.43 (?), with TMTpro tags manually added 515 as fixed modifications and used for quantitation The correction factors for purity of isotopic labels was set 516 according to vendor specification and minimum reporter precursor intensity fraction was set to 0.5. The 517 resulting protein groups were filtered for potential protein contaminants, protein groups only identified via 518 peptides decorated with modification or hits in the pseudo-reverse database used for FDR control. The 519 resulting intensities of each sample channel were normalised to the intensity of the 16th reference channel, 520 then median-centered and normalised according to the median-absolute deviation. Identified phosphopep-521 tides were similarly filtered, with the exception of filtering based on modified sites, and normalised using the 522 same strategy.

523
Differentially expressed phosphopeptides were called using the limma package (?) with a false discovery 524 rate of 0.05 on treatment minus control constrasts. Synergies were computed using a contrast fit of the 525 combination minus the sum of single treatments. Kinase substrate activity was implemented in R using the 526 ratio of the mean z-score as described in ? and computed for kinase-substrate sets from PhosphoSitePlus (?).

527
The normalised intensities and scripts used for the analysis can be found at https://itbgit.biologie.  (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted June 14, 2021. ; https://doi.org/10.1101/2021.06.14.448322 doi: bioRxiv preprint the IC50 data, the synergy data, the sequencing data and the proteomics data, and built the models. Falk Joern Toedling and Eric Blanc analysed the WES-Seq data. Joern Toedling and Clemens Messerschmidt 543 analysed the RNA-Seq data. Tommaso Mari measured and analysed the proteomics and phosphoproteomics 544 data. Anja Sieber performed the Western blot measurements. CC-BY-NC 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted June 14, 2021. ; https://doi.org/10.1101/2021.06.14.448322 doi: bioRxiv preprint