Bernardo-Faura et al. Genome Medicine (2021) 13:117 https://doi.org/10.1186/s13073-021-00925-8 RESEARCH Open Access Prediction of combination therapies based on topological modeling of the immune signaling network in multiple sclerosis Marti Bernardo-Faura1,2, Melanie Rinas3†, Jakob Wirbel1,3†, Inna Pertsovskaya4, Vicky Pliaka5, Dimitris E. Messinis6, Gemma Vila4, Theodore Sakellaropoulos5, Wolfgang Faigle7, Pernilla Stridh8, Janina R. Behrens9, Tomas Olsson8, Roland Martin7, Friedemann Paul9, Leonidas G. Alexopoulos5,6*, Pablo Villoslada4* and Julio Saez-Rodriguez1,3,10* Abstract Background: Multiple sclerosis (MS) is a major health problem, leading to a significant disability and patient suffering. Although chronic activation of the immune system is a hallmark of the disease, its pathogenesis is poorly understood, while current treatments only ameliorate the disease and may produce severe side effects. Methods: Here, we applied a network-based modeling approach based on phosphoproteomic data to uncover the differential activation in signaling wiring between healthy donors, untreated patients, and those under different treatments. Based in the patient-specific networks, we aimed to create a new approach to identify drug combinations that revert signaling to a healthy-like state. We performed ex vivo multiplexed phosphoproteomic assays upon perturbations with multiple drugs and ligands in primary immune cells from 169 subjects (MS patients, n=129 and matched healthy controls, n=40). Patients were either untreated or treated with fingolimod, natalizumab, interferon-β, glatiramer acetate, or the experimental therapy epigallocatechin gallate (EGCG). We generated for each donor a dynamic logic model by fitting a bespoke literature-derived network of MS-related pathways to the perturbation data. Last, we developed an approach based on network topology to identify deregulated interactions whose activity could be reverted to a “healthy-like” status by combination therapy. The experimental autoimmune encephalomyelitis (EAE) mouse model of MS was used to validate the prediction of combination therapies. * Correspondence: leo@mail.ntua.gr; pvillosl@clinic.cat; julio.saez@bioquant.uni-heidelberg.de † Melanie Rinas and Jakob Wirbel authors contributed equally to this work and are listed alphabetically. 5 School of Mechanical Engineering, National Technical University of Athens, Zografou, Greece 4 Institut d’ Investigacions Biomèdiques August Pi Sunyer (IDIBAPS), Barcelona, Spain 1 European Molecular Biology Laboratory, European Bioinformatics Institute, Cambridge, UK Full list of author information is available at the end of the article © The Author(s). 2021 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/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Bernardo-Faura et al. Genome Medicine (2021) 13:117 Page 2 of 16 Results: Analysis of the models uncovered features of healthy-, disease-, and drug-specific signaling networks. We predicted several combinations with approved MS drugs that could revert signaling to a healthy-like state. Specifically, TGF-β activated kinase 1 (TAK1) kinase, involved in Transforming growth factor β-1 proprotein (TGF-β), Toll-like receptor, B cell receptor, and response to inflammation pathways, was found to be highly deregulated and co-druggable with all MS drugs studied. One of these predicted combinations, fingolimod with a TAK1 inhibitor, was validated in an animal model of MS. Conclusions: Our approach based on donor-specific signaling networks enables prediction of targets for combination therapy for MS and other complex diseases. Keywords: Signaling networks, Pathways, Network modeling, Logic modeling, Kinases, Treatment, Personalized medicine, Combination therapy, Multiple sclerosis, Immunotherapy, Phosphoproteomics, xMAP assay Background combinations at the molecular level and allow us to pre- The signal transduction machinery is frequently affected by dict both efficacy and reduction of off-target effects [5, 11, perturbations induced by complex diseases. Hence, treat- 13]. ment of diseases such as cancer, cardiovascular, immuno- Multiple sclerosis (MS) is an autoimmune disease, in logical, or brain diseases is nowadays largely attempted by which the immune system is chronically activated and modulating different molecular cascades involved in the damages the central nervous system (CNS) [14]. The in- disease to stop its progression. As such, kinases involved in volvement of immune system deregulation in MS is signaling processes have evolved as primary targets for shown by altered phenotype and activity of blood lym- many diseases [1]. Further, there has been modest progress phocytes and monocytes [14], as well as by association with treatments based on single drugs. Combining several with genetic polymorphisms of immune genes [15]. At drugs targeting different pathways promises more effective present, there are fifteen Food and Drug Administration modulation of the pathogenic process [2, 3]. However, de- (FDA)-approved immunomodulatory drugs, and many velopment of combination therapies is hampered by the others are in late-stage clinical development. Most of often incomplete understanding on how their effect propa- these control the inflammatory activity in patients with gates through complex signaling networks, with crosstalk MS to a certain degree, though with known unwanted between the pathways influenced by each therapy [4, 5]. As effects at the signaling level. Furthermore, many unmet an additional level of complication, disease heterogeneity medical needs remain in the attempt of achieving con- hinders predicting how a specific combination therapy trol of the disease including more effective therapies, a could be translated into the clinic. Last, the combinatorial good safety profile, and neuroprotective or regenerative nature of such studies in terms of number of targets, drugs, treatments. Drug combinations are considered as a doses, and therapeutic regimens implies a large number of promising strategy to overcome some of these limita- experiments and associated costs, preventing a complete tions, and in cancer combination therapies, they are well analysis for all alternatives. As a result, the full potential of established [16]. However, predicting which patients combination therapies has not been fully developed yet. would benefit the most from a certain combination ther- Systems biology, and more specifically modeling of sig- apy remains as an unresolved challenge [17–20]. naling pathways applied to drug discovery, may provide a In this study, we present a systems medicine approach new path to approach this question [2, 6, 7]. Mechanistic aimed to (i) characterize the signaling pathways in pri- understanding at the network level offers integrated in- mary immune cells obtained from the blood of MS pa- sights about the cellular responses to environmental tients and healthy controls and (ii) predict new changes and drug effects, yielding a significant under- combination therapies based on the differences in signal- standing of the signaling cascades derived from decades of ing networks between treated MS patients and controls research in this field [8–10]. Mathematical modeling of (Fig. 1). To this end, we assembled a literature- and signaling networks has been used to unravel signaling database-based prior knowledge signaling network mechanisms and discover drug targets and disease media- (PKN) [21], which includes the pathways involved in im- tors such as cell surface receptors or intracellular mole- mune and MS signaling, as well as their crosstalk and cules by training those models to experimental in vitro known therapeutic targets. Based on their ability to measurements of key pathway components using inhibi- model large networks with a low number of parameters, tors and activators [4, 11, 12]. Modeling-based studies logic models have been used to unravel the network top- may be the key to characterize the effect of drug ology driving disease and response to therapy [22–25]. Bernardo-Faura et al. Genome Medicine (2021) 13:117 Page 3 of 16 Fig. 1 Topological modeling approach of signaling pathways for prediction of combination therapy. a Identification of subgroup networks. A model characterizing signaling activity (upstream kinase (circle) that regulates a response, e.g., innate immunity, survival (diamond)) in response to a stimulus (oval) was calculated for each donor based on the experimentally acquired dataset. Next, the donor-specific models are merged for all donors belonging to the same subgroup (left panel blue: healthy controls; middle panel orange: untreated MS; right panel green: treated MS). b Scoring subgroup interactions to find co-druggable network interactions. The score is calculated to identify interactions that differ from healthy- like signaling activity in spite of drug treatment (see the “Methods” section). c Topological prediction of drug combination. A topology-based graph search allows identifying secondary treatments that could target and revert signaling of co-druggable interactions to a healthy-like activity state Here, we established a Boolean logic model from the sig- lead to a personalized combination therapy. To identify naling network and trained it with measurements of kin- those interactions, we developed a score of co- ase de/phosphorylation as a proxy of signal propagation druggability of signaling interactions according to quan- upon perturbation with ligands and drugs in the periph- titative differences in network topology among healthy eral blood mononuclear cells (PBMCs) from MS patients controls and untreated and treated MS patients (Fig. 1 and healthy controls. We determined disease- and and Table 1). therapy-specific logic models that characterize the sig- Using this network-based approach, we predicted sev- naling networks for approved treatments (interferon- eral combination therapies, in particular of TAK1 with beta (IFNβ), glatiramer acetate (GA), natalizumab all five studied MS drugs. We validated a highly scoring (NTZ), fingolimod (FTY), and the experimental drug combination therapy in the animal model of MS experi- epigallocatechin gallate (EGCG)). We hypothesized that mental autoimmune encephalomyelitis (EAE). The mod- the signaling interactions that the drugs failed to revert eling approach shown here can be used for designing in the ex vivo assays to a healthy-like activity level may combination therapies for other complex diseases as well be candidates to be targeted by a second drug and hence as for developing personalized therapies. Bernardo-Faura et al. Genome Medicine (2021) 13:117 Page 4 of 16 Table 1 Overview of the identification of co-druggable interactions: co-druggability score examples and their interpretation Healthy MS Treated Score Treatment Affecting? Healthy-like signaling Interpretation after treatment 0 0 0 0 NO YES Healthy-like signaling: No combination therapy needed 0 0 1 -1 YES NO Un-healthy signaling. Single treatment is affecting: Co-druggable 0 1 0 1 YES YES Healthy-like signaling. Effective single treatment: No combination therapy needed 0 1 1 0 NO NO Un-healthy signaling. Single treatment not affecting: Co-druggable 1 0 0 0 NO NO Un-healthy signaling. Single treatment not affecting: Co-druggable 1 0 1 1 YES YES Healthy-like signaling. Effective single treatment: No combination therapy needed 1 1 0 -1 YES NO Un-healthy like signaling. Single treatment is affecting: Co-druggable 1 1 1 0 NO YES Healthy-like signaling. No combination therapy needed Co-druggable: interactions where treatment with the drug alone yielded signaling activity different to that of the healthy-like state. Columns 1–3: Signaling activity of a given interaction as assessed by modeling in healthy, untreated MS and treated MS. Column 4: Co-druggability score calculated based on differences between columns 1-3. Column 5: Difference present between treated and MS signaling. Column 6: Absence of a difference in interaction between healthy and treated signaling. Column 7: Combination of columns 4, 5 and 6 to identify interactions deregulated (i) by the disease and not reverted to healthy state by treatment or (ii) by an off-target signaling effect by the primary drug (see “Methods”) Methods Table 2 Demographic and clinical variables of MS patients and Subjects and clinical cohorts healthy controls (HC) We recruited 255 subjects including 195 patients with MS HC MS and 60 healthy controls in a multi-centric study in n=195 n=60 four MS centers (Hospital Clinic Barcelona – IDIBAPS Sex (M/F) 66/129 21/39 (n=69), Karolinska Institute (n=64), University of Zurich Age 43.1±11.3 39.9±8.5 (n=40), and Charité University (n=82)). Healthy donors, Disease duration (months) 104.9+93.2 -- untreated patients and patients treated with IFNB, GA, Age at onset 34.5+10.3 -- NTZ, and FTY were described before in an accompany- EDSS 2 (0-6.0) -- ing study [26], whereas the participants in a trial testing epigallocatechin-gallate (EGCG) at Charité University Disease subtype were described in [27]. Characteristics of the clinical co- CIS 24 -- horts: mean age: 43.1+11.3 years; disease duration: 8.7+ RRMS 129 -- 7.7 years; median Expanded Disability Status Scale SPMS 6 -- (EDSS): 2.0 (0–6.0); disease subtype: 24 clinically isolated PPMS 36 -- syndrome (CIS), 129 relapsing-remitting MS (RRMS), 6 Treatments secondary-progressive MS (SPMS), and 36 primary- progressive MS (PPMS); untreated: 93 and 60 healthy DMD controls (Table 2). In order to account for sex and age, FNb 36 -- healthy controls were matched to the RRMS subgroup, GA 18 -- i.e., the most prevalent subpopulation which is also the NTZ 22 -- MS subtype most frequently treated, a requirement to FTY 20 -- enable our goal to characterize signaling deregulation Experimental drug upon treatment. After quality control of the phospho- proteomic dataset, the number of subjects was reduced EGCG 6 -- to 169 (129 MS patients and 40 healthy controls), see Untreated Additional file 1: Supplementary methods. Untreated 93 -- M male, F female, EDSS Expanded Disability Status Scale, CIS clinically isolated Samples and processing syndrome, RRMS relapsing-remitting MS, SPMS secondary-progressive MS, PPMS primary-progressive MS, DMD disease-modifying drug, IFNb interferon A unified standard operating procedure for PBMC isola- beta, GA glatiramer acetate, NTZ natalizumab, FTY fingolimod, tion, stimulation, and lysis, as well as sample storing and EGCG epigallocatechin-gallate Bernardo-Faura et al. Genome Medicine (2021) 13:117 Page 5 of 16 shipping was developed along with a kit (plates) with Manual curation to prioritize interactions that were reagents and buffers that were produced in a single found in immune or human cells was used whenever facility (ProtAtOnce) and shipped to all participating possible. The global network featured (167 nodes and centers (see Additional file 1: Supplementary methods for 294 interactions). The identifiable components, i.e., details). The reagents were prepared from a single batch, those that led to unique solutions using the optimization and plates were prepared from a single batch for each algorithm, were assessed as described in [21], yielding a stimulus. Quality controls were carried out to ensure that so-called preprocessed network of 71 nodes and 168 in- the reagents remain stable for 3 months. teractions. To reduce model complexity to a degree that could be solved by model optimization, AND gates were xMAP assays added based on publicly available knowledge on protein XMAP assays were developed by ProtAtOnce (Athens, complexes at EBI’s complex portal [34] (https://www.ebi. Greece) and were standardized to minimize error. We ac.uk/intact/complex) (see Additional file 1: Supplemen- optimized assays from a list of 70 candidates (see Add- tary methods). Next, a personalized model for each donor itional file 1: Supplementary methods) and obtained a was generated using individual datasets as described in the final list of 17 phosphoproteins which display a good “Model optimization” section. signal to noise ratio to be measured in the in vitro assays: AKT1, CREB1, FAK1, GSK3A, HSPB1, IKBA, Model optimization JUN, MK03, MK12, MP2K1, PTN11, STAT1, STAT3, For each patient, 10 completed optimization runs were STAT5A, STAT6, TF65, and WNK1 (Additional file 2: performed with a genetic algorithm and assessed as Tables S1 and S2). We used a set of 20 stimuli, shown in Additional file 3: Figure S1 (see Additional which included pro-inflammatory or pro-oxidant sti- file 1: Supplementary methods). muli (Anti-CD3, concanavaline A (conA), IFNG, IL1A, IL6, LPS, NaCl, PolyIC, TNFα), immunomodulatory Network topology-based prediction of combination therapy stimuli (S1P, vitD3) neuroprotectants or anti-oxidants We calculated the co-druggability score as described in (BDNF, EGCG, INS, and BN201), disease-modifying Fig. 1 with S indicating the signaling activity status of drugs from MS (DMF, FTY, Teriflunomide, IFNβ1a each interaction: (Rebif®)), CNS-damaging oxidative stress H2O2, and a culture media as control (Additional file 2: Table S3). j SHealthy −S MS j − j SHealthy −S Treatment j Samples were collected at baseline (time 0) and after 5 were healthy refers to the healthy subgroup model, and 25 min. MS relates to the untreated MS subgroup model, and treatment denotes one of the five approved drug-treated Data normalization MS subgroup models. Therefore, interactions with a After signal reading, data was normalized. To allow logic negative co-druggability score indicated a treatment ef- modeling, data was normalized between 0 and 1, extend- fect that produced signaling activity more different from ing via stringent statistics the normalization strategy pre- that found in healthy donors than without treatment sented in [21] (see Additional file 1: Supplementary and were selected as co-druggable. A co-druggability methods). The maximum between 5 and 25 measure- score of 0 indicated that there was no effect due to drug ments was selected to allow capturing signal transduc- treatment. From those cases, the interactions in which tion including late effects of network motifs such as signaling activity was different between healthy donors negative feedback loops. and treated patients were also selected as co-druggable, i.e., those where the drug alone failed to revert signaling Model generation to a healthy state. A model describing core immune signaling was estab- lished based on Saez-Rodriguez et al [28]. Next, the j SHealthy −S MS j − j SHealthy −S Treatment j ≤ 0 model was extended to further immune and MS-related pathways such as interferon response, B and T cell re- Thus, we used the algorithm described in Fig. 1 to ceptor signaling, cellular survival and apoptosis, inflam- define the co-druggability of all interactions within each mation, lipid signaling, innate immunity, and multi-drug treatment group after filtering for active signaling. To response (MDR) genes from the state-of-the-art data- ensure that interactions with positive scores close to bases such as KEGG [29], ScienceSignaling [30], and zero were also captured by our method (interactions Wikipathways [31]. To allow the inclusion of MS drugs, with a similar signaling activity between drug treatment the drug targets were included from ChEMBL [32] and and MS untreated), we defined a lower quartile thresh- Uniprot [33] in pathways as explained above via the ref- old around zero and collapsed almost zero scores into erences detailed in Additional file 2: Tables S4 and S5. zero. Further, we required co-druggable interactions to Bernardo-Faura et al. Genome Medicine (2021) 13:117 Page 6 of 16 be significantly (p<0.05) different between healthy and Results drug signaling by using the upper quartile as the differ- Multiplexed phosphoproteomic analyses in PBMCs from ence threshold. Applying both filters together allowed us MS patients and controls to identify co-druggable interactions that were deregu- To characterize the signaling networks involved in MS, we lated by the disease as well as the unwanted signaling created a Prior Knowledge Network (PKN) of biochemical effect by the primary drug (Additional file 3: Figure S2). interactions (kinase phosphorylation) reported to be in- volved in immune signaling associated with MS based on Experimental autoimmune encephalomyelitis omics and functional studies, as well as those pathway in- Female C57BL/6 mice from Harlan (8–12 weeks old) teractions targeted by MS drugs [17]. Our network in- were immunized subcutaneously in both hind pads with cludes pathways such as interferon response, B and T cell 300 μg of a myelin oligodendrocyte glycoprotein receptor signaling, cellular survival and apoptosis, inflam- (MOG35-55, Spikem, Florence, Italy) emulsified with 50 mation, lipid signaling, innate immunity, and multi-drug μg of Mycobacterium tuberculosis (H37Ra strain; Difco, response (MDR) genes (Additional file 3: Figure S4). To Detroit, MI, US) in incomplete Freund’s adjuvant as de- achieve this, we searched in the state-of-the-art databases scribed previously [35]. Mice were injected intraperito- for interactions that were reported in highly specific assays neally with Pertussis toxin (500 ng; Sigma, US) at the and prioritized experiments with human and PBMC cells time of immunization and 2 days later. Animals were (see the “Methods” section). Further, targets of MS drugs weighed and inspected for clinical signs of disease on a were included in the network via their crosstalk with im- daily basis by an observer blind to the treatments. The mune pathways. The PKN featured 167 signaling compo- severity of EAE was assessed on the following scale: 0= nents (Additional file 2: Tables S5) and 294 interactions normal; 0.5= mild limp tail; 1= limp tail; 2= mild para- (Additional file 2: table S4). paresis of the hind limbs, unsteady gait; 3= moderate Subsequently, we developed a multiplex xMAP paraparesis, voluntary movements still possible; 4= para- phosphoprotein panel as a proxy of protein activation. plegia or tetraparesis; 5= moribund state; 6= death. The phosphosites in the panel were selected to Animals (10 mice per arm) were randomized to their maximize the coverage of the MS-specific network. The treatment once they have reached a clinical score > 1 stimuli were selected according to their regulation of the point. Clinical assessment was performed by a blinded pathways included in the PKN as detailed above and ref- evaluator to the treatment group. Comparison between erenced in Additional file 2: Tables S4 and S5. The groups was performed using the Mann-Whitney test panel, which combined previously and newly developed with p value cut-off at 0.05. assays, was then optimized to maximize accuracy, repro- ducibility, and network coverage (see Additional file 2: Statistical confirmation of patient subgroup models Table S1 and Additional file 1: Supplementary methods) To test whether the signaling models found for each [36]. After optimization, we selected a set of 17 phos- patient subgroup were reflected in the experimental phoproteins with adequate signal to noise ratio that were phospho-levels, we followed a three-stage strategy. First, used for the in vitro assays: AKT1, CREB1, FAK1, we identified stimulus-responding readouts in our experi- GSK3A, HSPB1, IKBA, JUN, MK12, MK03, MP2K1, mental dataset by assessing for each stimulus-readout PTN11, STAT1, STAT3, STAT5, STAT6, TF65, and combination if the fold changes (before non-linear WNK1 (Additional file 2: Table S2). PBMCs were cul- normalization) of a given readout deviated significantly tured in the presence of different sets of stimuli such as from 0 for all donors within a single patient subgroup lectins, endotoxins, immunostimulants, cytokines, and using a one-sample Wilcoxon test. Multiple testing cor- drugs (Additional file 2: Table S3). rection was performed using the Benjamini-Hochberg To estimate the representativeness of our selected strategy, yielding a p value for each stimulus-readout com- phospho-signals, we calculated closeness centrality, a bination within a group. Secondly, we determined which metric which measures the efficiency of a node in stimulus-readout pairs were either in accordance with the spreading information through a graph, for all nodes in model of signaling or not, by confirming if a phospho- the immune- and MS-specific signaling network pre- readout could be reached from the respective stimulus in sented here. No signal was found to be central to the the signaling model for the different subgroups whole network, and the nodes we selected and measured (Additional file 3: Figure S3A-E). Finally, using Fisher’s as signals were found to be the most central nodes in exact test, we calculated whether the stimulus-readout the B cell, T cell, inflammation, IFNb, cellular survival, pairs found to be significant (using p<0.05 as cutoff) in the apoptosis, and innate immunity pathways (Additional first stage were enriched for those supported by the model file 3: Figure S5), thereby supporting that they were (determined in the second stage) compared to those not indeed representative of the pathways we sought to supported (Additional file 3: Figure S3F). study. Bernardo-Faura et al. Genome Medicine (2021) 13:117 Page 7 of 16 We extracted PBMCs from 195 MS patients (mean samples (see Additional file 1: Supplementary methods), age: 43.1+11.3 years; disease duration: 8.7+7.7 years; me- which led to the removal of 2028 data points and 11 pa- dian Expanded Disability Status Scale (EDSS): 2.0 (0– tients to generate the final dataset for a total of 169 sub- 6.0); subtype: 24 CIS, 129 RRMS, 6 Secondary- jects. To reduce the complexity of the model to those Progressive MS (SPMS), and 36 PPMS; untreated: 93) pathways that can be determined based on the experimen- and 60 healthy controls (Table 2). From those 255 do- tal coverage, identifiability analysis was performed as de- nors, 180 were selected and PBMCs were analyzed be- scribed in [21]. After identifiability analysis, the PKN fore stimulation and at 5 and 25 min after stimulation. (shown in Additional file 3: Figure S4) was reduced to 71 This yielded a dataset consisting of three measurements proteins and immune modulators and 168 interactions. for 17 phosphoproteins upon 20 stimuli, to a total of 183,600 experimental measurements. For subsequent Logic modeling enables characterization of signaling modeling, the two-time points after baseline (5 and 25 networks in healthy donors and MS patients min) were collapsed into a single activation signature to Next, we sought to train the generic PKN to these phos- define whether a given phosphoprotein was activated or phoproteomic measurements in a donor-specific manner inhibited (see details in the “Methods” section). Finally, to identify each donor’s active pathways. To train the we applied stringent quality control analyses of the data- same Boolean logic model to the data of each of the 169 set consisting of positive and negative controls, reproduci- subjects, the phosphoproteomic datasets were trans- bility tests, and tests for artifacts in the distribution of the formed (Fig. 2a). To enable comparison of the data to Fig. 2 Phosphoproteomic measurement and normalization pipeline. a xMAP mean fluorescence intensity (MFI) log values of the 17 analyzed phosphoproteins, b fold change distribution, c non-linearly normalized values (see the “Methods” section). Orange measurements a–c: values of the same patient to allow visualization of the changes across data transformation. d Percentage of patients, for which each phosphoprotein was classified as phosphorylated, dephosphorylated, or non-significant after statistical testing Bernardo-Faura et al. Genome Medicine (2021) 13:117 Page 8 of 16 the binary output of the Boolean model, we established a than MS untreated patients (median Jaccard distance = stringent normalization pipeline combining (i) a non- 0.644) to the mean healthy signaling network. We then linear transformation to normalize the measurements to calculated pairwise distances among all donors, and used continuous values between 0 and 1 while penalizing the them as a background similarity for comparison. We outliers [21] with (ii) a statistical filter that removed found the dissimilarity among healthy donors to be sig- data-points belonging to proteins which, upon perturb- nificantly lower than that of the background (Wilcoxon ation, were not found to be significantly phosphorylated test’s p value = 0.0215). The dissimilarity in untreated or dephosphorylated. Following this strategy to identify patients was also lower than the background’s (p value = two significant states from continuous phosphorylation 0.0318). In addition, drug treatment exhibited a strong measurements featuring a three-state scenario of phos- effect on signaling, which seemed to be homogenized phorylation, dephosphorylation, or unchanged levels, within each group with dissimilarities lower than that of this procedure enables us to train a Boolean (binary) healthy and untreated donors (median Jaccard distance model (Additional file 1: Supplementary methods, “Data to each group’s mean network: EGCG = 0.512, untreated normalization for logic modeling” and Additional file 3: MS = 0.641; FTY = 0.595, untreated MS = 0.646; GA = Figure S6). Figure 2b shows the log fold change of each 0.632, untreated MS = 0.665; NTZ = 0.615, untreated phospho-profile with respect to the control, Fig. 2c the MS = 0.641; IFN = 0.599, untreated MS = 0.615). The dataset after the non-linear normalization, and Fig. 2d differences were found to be significant in three of the displays the proteins found to be significantly dephos- treatments (Wicoxon’s test p value EGCG = 0.00676, phorylated or phosphorylated. Using this normalized FTY = 0.000418, IFNb = 0.0462), further supporting the dataset, we aimed to identify the active specific signaling strong effect of drug treatment on signaling. Altogether, network for each donor. these results supported that merging models by sub- To this end, we fitted a general logic model derived groups of donors yields biologically meaningful signaling from the PKN [21] to each donor dataset, which resulted networks for each group. in a compendium of 169 models. Model optimization Next, we sought to characterize the specific network ar- was performed with CellNOpt, a tool that selects the chitectures of MS, as well as those of MS treatments. The logic model that best matches the data while penalizing differentially activated interactions between healthy donors model size [37]. To increase the robustness of each and untreated MS patients uncovered signaling pathways model, optimization was performed 10 times, and for deregulated in MS (Additional file 1: Supplementary each donor, we selected the best solution as well as all methods and Fig. 3c). Our analysis revealed the activation models with a data-to-simulation mismatch within a of several pathways in PBMCs after stimulation both in pa- given relative tolerance reflecting that, due to lack of tients and controls, ranging from cell survival and prolifera- identifiability and technical noise, different models are tion to TCR, innate immunity, and pro-inflammatory feasible [21] (see the “Methods” for details). We then response pathways (e.g., TCR - CSK - LCK, JAK1 - STAT1, built a final model for each patient, defined by the me- TLR4/IL1R1 - MYD88 - IRAK4 - TRAF6 - TAK1, TNFα - dian value of each interaction across all solutions within TRADD - TRAF2 - MAP3K1, SGK - MK03, or RS6 - SGK, the relative tolerance (or borders of tolerance) of the RAF1 - MP2K1). Further, increased activity was found in best solution (Fig. 3a). To assess the validity and robust- pathways such as INSR - PI3K - PIP3 - AKT1 and JAK1 - ness of our modeling approach, we confirmed the ab- STAT3 in MS patients, while BDNF - NTRK1 - GRB2 was sence of bias due to treatment, center, disease subtype, found to be activated in healthy donors. medical condition, and technical aspects such as model These findings indicated that our method was able to size and merging strategy (Additional file 3: Figure S1). identify previously described pathways in the setting of To study the effect of MS and MS treatment on signal ex vivo analysis of human PBMCs [17], as examined in transduction, we first merged the individual donor detail for the individual pathways in the discussion. Un- models within each group by calculating the mean activ- treated MS patients, when compared to healthy controls, ity over each reaction, which yielded a single network showed an increase in the activation of the NFkβ path- per group. Next, we calculated the distance between way (TAK1-IKKB), the activation of the cell prosurvival each patient model and the corresponding group mean PI3K pathway (SLP76-AKT1), or the activation of inter- network using the Jaccard distance, a metric used to as- feron/cytokines pathways (JAK1-STAT3). Moreover, pa- sess the number of interactions that are dissimilar in tients showed a decrease, compared to healthy controls, two networks [38]. Finally, as a comparison, we calcu- in the activation of TCR/IL-2 pathway (LCK-STAT1) as lated for each group the distances between the mean well as in the effect of trophic factors signaling on ubi- group network and the networks of all untreated MS pa- quitination system (EGFR-CBL). These results tients (Fig. 3b). We found that healthy donors were characterize, at the mechanistic and quantitative level, slightly more similar (median Jaccard distance = 0.628) the differential activation of the immune pathways in Bernardo-Faura et al. Genome Medicine (2021) 13:117 Page 9 of 16 Fig. 3 Logic modeling identifies donor-specific signaling networks and reveals MS-specific signaling pathways. a Signaling network found by modeling for each donor, visualized as a heatmap. Rows: Single donor network. Columns: Signaling activity determined for each interaction by calibrating the PKN shown in Additional file 3: Figure S4 after removing the unidentifiable interactions using the phosphoproteomics dataset of each donor. b After networks were merged by subgroup, the Jaccard distance was used to assess similarity from all donors within each group (selected donors in group legend) to their mean subgroup network (network in X axis) and compare it to the similarity from MS patients to the same group network. Healthy donors (blue) were more similar to the mean healthy network than untreated MS patients (orange). In turn, the distance from both groups of donors to that of the combined signaling activity in all donors (grey) was statistically significant. Distance from treated donors (green) to their mean subgroup network was largely reduced when compared to distance from untreated donors to the treatment’s network, suggesting a strong effect of treatment homogenizing within group signaling. c Differentially activated pathways (see Additional file 1: Supplementary methods) between healthy controls (HC) and untreated MS patients (MS). The models previously calculated for each donor were merged to reveal the common active pathways for controls (blue), untreated MS patients (orange), and both (brown). Gray: Inactive interactions from the MS, immune- and treatment-related network (Additional file 3: Figure S4) Bernardo-Faura et al. Genome Medicine (2021) 13:117 Page 10 of 16 MS and were supported by the statistical analysis of the in spite of its known efficacy in MS, yielded a signaling phosphorylation of individual kinases in an accompany- activity that differed from the healthy controls as identi- ing study [26]). In particular, the accompanying study fied by our topological modeling approach. Therefore, found MP2K1, STAT1, STAT3, TF65, and HSPB1 to be we aimed to identify which kinase interactions within differentially phosphorylated (p<0.05 after correction for the network needed its signaling activity to be reverted multiple hypothesis testing) in MS patients relative to the to the healthy-like state when treated with the ongoing controls, confirming the differential activation of cell therapy. We hypothesized that co-druggable interac- survival and proliferation (MAPK), and pro-inflammatory tions, i.e., those that a given MS therapy failed to revert (STAT) pathways in immune cells. to a healthy-like activity level, should be the model inter- actions with a signaling value more distant between Quantitative differences found in signaling pathways healthy control and MS drug models than between under therapy healthy control and untreated MS models (Fig. 1, Table We then quantified differences in signaling between pa- 1 and see the “Methods” section). This yielded a list of tients treated with MS-specific therapies, using the same interactions with their corresponding score quantifying procedure as for healthy donors versus untreated MS pa- co-druggability potential (Additional file 2: Table S6 and tients. We focused on the strongest signals, i.e., those in Additional file 3: Figure S2). To identify the co- the upper quartile of the group mean (Additional file 3: druggable interactions, we selected those with a non- Figure S2). This yielded a signaling network of active positive score as described above and filtered for add- reactions characterizing each subgroup, uncovering the itional conditions (see the “Methods” section), thereby effect of each MS-treatment on signaling at the mechan- identifying both the interactions whose signaling activity istic level (Additional file 3: Figure S3). was different from healthy because of deregulation by Interestingly, the activation of TAK1, a key component the disease as well as because of off-target primary drug involved in TGF-β, Toll-like receptor, B cell receptor, effect (Fig. 4a). The last step of our approach was to pre- and NFkβ pathways, as well as response to inflammation dict the stimulus that would revert to healthy-like signal- was identified in all five networks inferred under treat- ing activity of those interactions identified as co- ment (Additional file 3: Figure S3 A-E). Additionally, its druggable. Therefore, the co-druggable interactions were activity was persistently over-activated in terms of quan- mapped onto the signaling network assessed for each titative signaling activity by the currently approved MS drug subgroup. As an example, Fig. 4b shows the co- drugs (Additional file 3: Figure S2). In the case of IFNβ- druggable interactions under FTY treatment mapped and GA-treated patients, the network revealed the acti- onto the FTY network. Mapping the co-druggability for vation of STAT1 and STAT3 by JAK1, whereas in the each interaction allowed us to predict combination ther- EGCG-, FTY-, and NTZ-specific networks, the activa- apies based on each treatment’s signaling network top- tion of the JAK1-STAT pathway is regulated via STAT3, ology. Finally, we employed a graph search approach to STAT6, and STAT1, respectively. identify the in vitro stimulus used in our study that acti- To validate that the therapy-specific pathways found vated interactions found to be co-druggable with each were consistent with the experimental measurements of drug. In other words, we found the experimental read- MS signatures, we quantified the degree to which the outs that could be measured and reached from an original phospho-measurements supported the signaling in vitro stimulus via an interaction found to be co- pathways predicted by modeling for each patient sub- druggable in vivo (Additional file 2: Table S7). It is im- group. The proteins differentially phosphorylated across portant to note that our approach enabled the prediction all stimuli combinations were overrepresented in path- of combinations of approved MS therapies which yielded ways found by model fitting (IFNβ p=1.2E−05; GA p= a signaling activity distant from healthy controls with 1.4E−06; FTY p=0.0004; NTZ p=0.006; Fisher’s exact drugs that stimulated or inhibited signaling. We statisti- test) except for EGCG, which was limited by the small cally confirmed that the signaling models found for each sample size (Additional file 3: Figure S3F and methods). patient subgroup were reflected in the experimental phospho-levels (see statistical confirmation in the “Methods” Network topology-based prediction of targeted section and Additional file 3: Figure S3). combination therapy Our main goal was to use the subgroup networks found Validation of a predicted combination therapy in the to predict novel combination therapies. To this aim, we animal model of MS defined a therapeutic goal to revert the signaling net- Next, we validated in vivo one of our predictions. work of patients with MS to a healthy-like activity. We Among all predictions (Additional file 2: Table S7), we introduced a restriction to our method: the combina- chose to validate the combination of FTY with a TAK1 tions we sought included one approved MS drug which, inhibitor based on (i) the striking signaling homogeneity Bernardo-Faura et al. Genome Medicine (2021) 13:117 Page 11 of 16 Fig. 4 Combination therapies predicted and in vivo validation. a All predicted co-druggable interactions of the MS drugs models. Based on the subgroup models, the co-druggability of all 168 network interactions (X axis) was assessed using the co-druggability score, and those identified as co-druggable (see Fig. 1, Table 1 and main text) are shown. For each interaction (X-axis) the number of drugs (Y-axis) is shown, in which it was found to be co-druggable using the co-druggability criteria. b FTY network co-druggability: the case of FTY network co-druggability is shown as an example (red line: interactions predicted to be co-druggable). c In vivo validation of the combination FTY+TAK1-inhibitor in the EAE model. The graph shows the mean and the standard error of the clinical score for each group (n=7). Animals started treatment after disease onset (clinical score >1.0) and were randomized to each treatment and rated in a blinded manner. Stars show days significantly different from the same day with placebo across FTY models (as quantified in Fig. 3b), (ii) the EAE model, whereby the next lower dose without effi- finding that TAK1 modulation of HSPB1 via MK12 is cacy below the efficacious one was selected. We found largely deregulated (Additional file 3: Figure S2, panel that the combination FTY+TAK1-inhibitor ameliorated FTY, interaction MK12 - HSPB1), and (iii) the fact that the clinical course of EAE compared to placebo, TAK1- TAK1 is active and co-druggable in all 5 networks under inhibitor, or FTY alone (Mann-Whitney test; p<0.05) treatment with in vivo drugs, specifically TAK1 - (Fig. 4c). Therefore, our method was able to identify MAP2K4 for ECGC, IL1R1 - TAK1 for FTY, IL1R1 - combinations of drugs targeting different pathways that TAK1 for GA, LPS - TAK1 for IFNβ, and TAK1 - MK12 achieved higher efficacy than single therapy. for NTZ (Additional file 2: Table S6). To obtain an in vivo proof-of-concept of the efficacy Discussion of the combination therapies proposed, we sought to as- In this study, we have built donor-specific dynamic logic sess whether the combination therapies improved the models and developed a network-based approach to (i) clinical course of the animal model of MS (EAE in the characterize signaling deregulation in MS and several of C57BL6 mice immunized with the MOG35-55 peptide). its treatments and (ii) predict new targets for combination Mice were randomized after disease onset (after they therapy. Although highly effective therapies have already reached a clinical score >1.0), either to placebo (saline), been developed for MS, improving their efficacy further, each drug alone, or the combination of FTY with the particularly for progressive MS, remains an important un- TAK1 inhibitor (5Z-7-oxozeaenol). Doses were selected met medical need [41]. Higher efficacy is desirable for from previous dose-efficacy studies [39] and [40] in the progressive MS, in which currently approved drugs Bernardo-Faura et al. Genome Medicine (2021) 13:117 Page 12 of 16 (Ocrelizumab and Siponimod) show approximately 20% supported by multiple previous studies identifying deregu- amelioration in primary and secondary progressive MS re- lation of pathways in MS [17] and shows the activation of spectively compared to placebo. Even in the case of pro-inflammatory pathways and the inhibition of path- relapsing-remitting MS, not all patients reach complete ways related with immune tolerance. control of the disease (no evidence of disease activity Furthermore, the phosphoproteomic dataset designed (NEDA)). Recent studies have found that RRMS patients for and used in this study was also used in an accom- accumulate disability independently of relapses even if panying work [26] that performs a statistical comparison treated with low or high efficacy therapies as shown in the of phosphorylation of single kinases and their association Ocrelizumab trials [42]. Whether the combinations that to genetic susceptibility. The accompanying study found we identified with our approach would improve such multiple differentially phosphorylated and associated outcome will require to be tested in clinical trials. kinases, which are crucial components of the pathways In addition to predicting combination therapies, the found deregulated here. approach presented here allows insights into signaling The combination therapies predicted with our topo- deregulation in MS by comparing signaling networks in logical modeling approach may contribute to reducing untreated MS patients with those of healthy controls. The the treatment’s unwanted effects in two ways. First, our comparison demonstrates enhanced pro-survival effects of method enables the characterization of drug-specific ef- the trophic factor signaling pathway (AKT1) and modula- fects on signaling. For example, the models revealed the tion of the interferon pathway (JAK1, STAT3). In line with activation of JAK1-STAT1 pathway [60, 61] in IFNβ- the differences found here, aberrant STAT phosphoryl- treated patients, the activation of AKT, PLC [62, 63] in ation signaling in peripheral blood mononuclear cells patients treated with FTY, or the activation of MAPK, from MS patients has been reported [43]. The NFkβ path- and NFkβ (via TAK1) [64–66] in NTZ-treated patients. way has been reported to be overactivated in PBMCs from By identifying models of MS drugs, our method found patients with MS [44, 45] as well as to contribute to the the pathways activated downstream of drug targets in genetic susceptibility of the disease [46, 47]. Our analysis our signaling network. Further, it identified interactions supports the systemic pro-inflammatory state of PBMCs with other disease-associated pathways, which may not in MS. The trophic factor pathway involving SLP76 and directly be involved in the signaling targeted by such a AKT was found active in patients under treatment with drug. Therefore, our modeling approach can be used to GA, NTZ, and IFNβ and has also been associated with provide a wider context of the effects of a given drug on MS [17] and MS susceptibility via CD6 gene [48, 49]. The the functions of the target cells. Specifically, our co- involvement may reflect a pro-survival signaling state of T druggability score identifies not only those interactions and B cells in the context of the pro-inflammatory micro- that are unhealthy-like due to MS deregulation but also environment. Finally, we observed overactivation of the those that change due to first line treatment, whose sig- cytokine/interferon pathway (JAK1), which has previously naling can be reverted to healthy-like by a combination been reported in MS [17]. Moreover, STAT3 has been drug. Second, rational identification of drug combina- confirmed as susceptibility gene for the disease [50], and tions may allow us to use lower doses of approved drugs, its activation is impaired in response to IL-10 in MS thereby reducing the dose-related adverse events of each patients [51], suggesting a defective response of regulatory single drug, improving safety and tolerability while main- Tr1 cells. Regarding the interactions that were decreased taining efficacy. As we show in this study, the latter in MS patients compared to controls, PBMCs from aspect can be predicted from in vitro measurements and patients with MS showed lowered inhibition of STAT5 by validated in animal models, as we did with the combin- LCK suggesting impairment of the regulation of T/B cell ation of FTY and a TAK1 inhibitor. In the specific case signaling and IL-2 trophic effects [52] or cytotoxicity [53], of FTY, its most common side effects are dose- as well as the regulation of the ubiquitination system dependent and include lymphopenia, transient bradycar- modulated by CBL-B [54]. LCK is modulated by EVI5, dia, increased rate of latent infections with HSV-1 and and influences STAT5 in our analysis, both being suscep- VZV, liver enzyme elevation, and a few others. Fingoli- tibility genes for MS [50]. The second inhibited interaction mod acts by trapping lymphocytes in secondary lymph- involved the MS susceptibility gene CBL-B [55], which oid organs by functional antagonism of S1P receptors. regulates TCR and co-stimulatory signals and regulates However, its effects on cell signaling are poorly under- immune tolerance through its ubiquitin E3-ligase activity. stood. Some of its adverse events may be prevented by CBL-B expression is reduced in CD4 cells from MS reverting signaling to healthy-like in the TGF-β, Toll- patients and alters the signaling of the type I interferon like receptor, B cell receptor, NFkβ, or proinflammatory pathway [56, 57]. CBL-B is activated by EGFR, leading to pathways of which TAK1 is a member. Further, in our inhibition of several pathways by ubiquitination, including validation, we made use of suboptimal doses and showed that of EGFR itself [58, 59]. In summary, our results are that the combination is effective in EAE mice as it is Bernardo-Faura et al. Genome Medicine (2021) 13:117 Page 13 of 16 with Fingolimod at full dose, although using significantly our models. Third, the Boolean logic approach does not de- smaller doses. Hence, safety can be reasonably expected scribe processes biochemically and models processes as to be improved. While these considerations regarding in- binary, hence missing subtle aspects of signal transduction creased safety by drug combinations are reasonable, [10]. One of the main goals of our study was to develop a demonstrating them is likely more difficult and will re- modeling approach that can guide the rational development quire careful testing in patients. For instance, to date of combination therapies. Our signaling topology-based ap- TAK1-inhibitors have not been and are not in clinical proach has the advantage of allowing prediction of combi- testing according to ClinicalTrials.gov for MS. However, nations between drugs currently used in patients with MS TAK-1 inhibitor Takinib has been shown to broaden the and compounds that can stimulate signaling where those therapeutic potential of anti-TNF approaches in cancer drugs yield a signaling activity distant from the healthy and autoimmune disease [67], and we anticipate that it state. Many strategies have been developed to predict drug is going to be developed for clinical testing. We hope combinations [71]. Others have used phosphorylation data that these insights can be used to reduce unwanted side upon perturbation using statistical approaches [13], data- effects and tailor treatment to patients. driven network inference [5, 11], or a combination of Our study has several limitations. First, the analysis of sig- mechanistic and Bayesian network modeling [72]. Due to naling signatures is based on phosphoproteomics of mixed their simplicity, logic networks can model large networks immune cells, i.e., PBMCs. While this allows us to measure and provide a useful framework to study drug combinations key pathogenic MS features with a focus on immune mech- [23–25]. We used this logic framework to analyze a large anisms such as T cell and B cell activation, cell adhesion dataset encompassing 183,600 data points derived from a and migration, proinflammatory differentiation, and others newly recruited cohort of 169 donors. Our approach is a in the multiple cell types that have been related to MS compromise between data availability, technical feasibility, pathogenesis [68], it yields a signal that is the average of the and computational burden that sacrifices details to be able signals in cell subtypes and hence can mask cell-type- to capture a broad portion of the signaling machinery. specific responses. To avoid a signal averaged across Fourth, in this study, we have explored only some of the PBMCs, one could filter specific immune cell subtypes such currently approved drugs for MS patients. And fifth, as CD4+, CD8+, B cells, and monocytes. Subsequently ana- we have averaged individual patients’ networks to ob- lyzing signaling abnormalities while accounting for differen- tain the subgroup (drug) networks, which leads to tial contribution to disease susceptibility or response to losing individual variability, probably related to differ- therapy across cell subtypes may reveal new therapeutic tar- ential genetic susceptibility, or immune system gets. Alternatively, technologies to collect information at activation. the single-cell level such as mass cytometry [69] would In summary, we have built donor-specific dynamic allow us to separate cells by type, but the cost is consider- logic models and developed a network-based approach ably higher, limiting the number of subjects. A more afford- to (i) characterize MS and treatment deregulation of sig- able alternative would be flow cytometry, which on the naling and (ii) predict new targets for combination ther- other hand can measure a lower number of analytes. Sec- apy. This approach can be applied to other diseases. ond, our current coverage using validated xMAP assays was restricted to 17 proteins. Although designed to maximize Conclusions the coverage of immune- and MS-related pathways and Here, we applied a modeling approach to characterize confirmed as representative by topological analysis, these signaling activity both at the donor-specific and at the constitute a small subset of the signaling molecules and subgroup level using phosphoproteomic data from pathways that may be activated in immune cells. Some primary immune cells of MS patients and healthy donors. technologies such as mass spectrometry allow studying At the subgroup level, we identified models of the signal- multiple phosphosites and at whole proteomic scale. How- ing network for healthy controls, untreated, and drug- ever, the specific proteins measured are not chosen and treated MS patients. Our approach allowed us to cost per sample is much larger. Therefore, the number of characterize signaling deregulation in this complex and patients, timepoints, and replicates would have to be re- heterogeneous disease. Further, based on the data-driven duced given our budget, making our study unfeasible. In models, we developed a network-based method to predict those cases where upscaling the number of readouts is pos- novel combination therapies for the treatment of MS. We sible, robust quantitative phosphoproteomic assays covering hypothesized that the interactions in the subgroup models thousands of phosphosites in hundreds of samples can be in response to single drugs that were not reverted to a used for logic modeling [70]. This will allow us to signifi- healthy-like activity state can be co-drugged. Hence, we cantly expand the scope of our models. Furthermore, other developed a co-druggability score identifying those inter- modifications such as ubiquitination can play an important actions. Finally, we used that score on a newly developed role in signaling. Such data, if available, can be included in strategy to predict combination therapies based on Bernardo-Faura et al. Genome Medicine (2021) 13:117 Page 14 of 16 network topology. Thus, our algorithm identifies interac- Funding tions within signaling networks that should be targeted to This work was supported by the European Union 7FP-Programme (CombiMS, grant No 305397) and the European Sys4MS project (Horizon2020: Eraco- restore the network to a healthy-like state. sysmed: ID-43). RM and the Section for Neuroimmunology and MS Research, As validation, we tested in vivo the combination of University Zurich (UZH), have been supported by a European Research Coun- one approved drug for MS, FTY, with an inhibitor of cil Advanced Grant (ERC 340733-HLA-DR15 in MS), by a Clinical Research Pri- ority Project - disease heterogeneity of MS (CRPPMS; UZH), by the Swiss TAK1, a key component in the predicted signaling National Science Foundation and the Swiss MS Society. Open Access funding networks of all five studied in vivo drugs, which was enabled and organized by Projekt DEAL. found highly deregulated in FTY. This combination Availability of data and materials largely ameliorated the course of the animal model of The data from clinical studies has been published by Kotelnikova et al. in MS with significantly higher efficacy than each treat- [26], while the data from the EGCG clinical trial has been published by ment alone, providing in vivo proof-of-concept for Bellmann-Strobl et al. in [27]. All data generated in this study are available on Github (http://github.com/saezlab/combiMS) [73], including all data transfor- our approach. mations from xMAP raw measurements to modelling-transformed as shown In conclusion, the approach developed here can be in Fig. 2. All code and mathematical models can be found under the same applied to other diseases with poorly understood patho- Github repository. genesis and treatments which may produce severe side Declarations effects or only partially ameliorate the disease. We ex- pect our approach will contribute to a better under- Ethics approval and consent to participate standing of disease- and therapy-related signaling Patients were invited to participate and included in the study by their neurologist in charge of medical care and included after signing the deregulation and aid in the development of combination informed consent. All appropriate IRB approvals to use the data from the therapies to restore healthy signaling. clinical studies and the EGCG trial were obtained from these IRB committees with the following committee reference numbers. University of Zurich: Approved by the Cantonal Ethics Board of the Canton of Zurich, EC-Nr. 2013- Abbreviations 0001. Karolinska Institutet: Approved by Regionala etikprövningsnämnden i MS: Multiple sclerosis; EGCG: Epigallocatechin gallate; TAK1: TGF-β activated Stockholm: Stop MS II, EC-Nr. 2009/2107-31/2. IDIBAPS: Approved by Comite kinase 1; TGF-β: Transforming growth factor beta-1 proprotein; CNS: Central de Etica Medica Hospital Clinic Barcelona, EC-Nr: 2012-7662. Charité Univer- nervous system; FDA: Food and Drug Administration; PKN: Prior Knowledge sity Medicine Berlin: Approved by Ethikkommission des Landes Berlin am Signaling Network; PBMC: Peripheral blood mononuclear cell; Landesamt für Gesundheit und Soziales, EC-Nr. 2006-006323-39. Readers can IFNβ: Interferon-beta; GA: Glatiramer acetate; NTZ: Natalizumab; gain access to the EGCG trial data as described in the original publication FTY: Fingolimod; EAE: Experimental autoimmune encephalomyelitis; [27] by contacting friedemann.paul@charite.de, as far as permitted according RRMS: Relapsing-remitting MS; SPMS: Secondary-progressive MS; to data protection requirements and consent provided by the participants. PPMS: Primary-progressive MS; CIS: Clinically isolated syndrome; The research conformed to the principles of the Helsinki Declaration. EDSS: Expanded Disability Status Scale All animal handling was carried out in accordance with the European Council Directive (2010/63/EU) and the Spanish regulations for the procurement and care of experimental animals (1201 RD/2005, October 10). Supplementary Information All the study protocols were approved by the Ethical Committee on Animal The online version contains supplementary material available at https://doi. Research of the University of Barcelona. org/10.1186/s13073-021-00925-8. Consent for publication Not applicable. Additional file 1. Supplementary methods for phosphoproteomics dataset generation, data processing, statistical analyses, mathematical modeling and combination therapy prediction. Competing interests DEM is an employee of ProtATonce; TO has received honoraria for lectures Additional file 2: Supplementary Tables S1 to S7. and/or advisory boards as well as unrestricted multiple sclerosis research Additional file 3: Supplementary figures S1 to S6. grants from Allmiral, Astrazeneca, Biogen, Genzyme, Merck, and Novartis; RM reports grants and personal compensation from Biogen, personal fees from Genzyme Sanofi Aventis, grants and personal fees from Novartis, and Acknowledgements personal fees from Merck Serono, Roche, Neuway, and CellProtect, outside The authors thank Federica Eduati for insightful discussion and valuable the submitted work; FP has received research grants and personal advice all along the way. We also would like to thank Thomas Cokelaer for compensation for activities with Alexion, Bayer, Chugai, Novartis, Merck, Teva, crucial technical help. Sanofi, Genzyme, Biogen, and MedImmune; LGA is the founder and hold stocks at ProtATonce; PV holds stocks and has received consultancy fees from Bionure Farma SL, Spiral Therapeutics Inc., Spire Bioventures Inc., Attune Authors’ contributions Neurosciences Inc., QMenta Inc., and Health Engineering SL; JSR declares Developed the modeling and combination therapy prediction approach: funding from GSK and Sanofi and consultant fees from Travere Therapeutics. MBF. Conceived and designed the experiments: MBF, IP, VP, DEM, LGA, JSR, The other authors declare that they have no competing interests. and PV. Developed new phospho-assays: VP, DEM, TS, and LGA. Designed pa- tient cohort and recruited patients: TO, RM, FP, and PV. Performed the experi- Author details 1 ments: VP, DEM, GV, WF, PS, and JB. Contributed reagents/materials/analysis European Molecular Biology Laboratory, European Bioinformatics Institute, tools: MBF, TO, RM, FP, LGA, JSR, and PV. Analyzed the data: MBF, JW, TS, and Cambridge, UK. 2Centre for Research in Agricultural Genomics (CRAG), JSR. Tested algorithm, performed statistical corroboration and additional ana- CSIC-IRTA-UAB-UB, Campus UAB, Bellaterra, Barcelona, Spain. 3Joint Research lyses, and modified figures accordingly: JW and MR. Wrote the manuscript: Center for Computational Biomedicine (JRC-COMBINE), Faculty of Medicine, MBF. Edited the manuscript: JSR and PV. Supervised the clinical study and ex- RWTH-Aachen University, Aachen, Germany. 4Institut d’ Investigacions perimental validation: PV. Supervised the modeling, combination therapy Biomèdiques August Pi Sunyer (IDIBAPS), Barcelona, Spain. 5School of prediction, and overall project: JSR. The authors read and approved the final Mechanical Engineering, National Technical University of Athens, Zografou, manuscript. Greece. 6ProtATonce Ltd., Athens, Greece. 7University of Zurich, Zurich, Bernardo-Faura et al. Genome Medicine (2021) 13:117 Page 15 of 16 Switzerland. 8Department of Clinical Neuroscience, Karolinska Institutet, 20. Milo R, Panitch H. Combination therapy in multiple sclerosis. J Stockholm, Sweden. 9NeuroCure Clinical Research Center and Department of Neuroimmunol. 2011;231(1-2):23–31. https://doi.org/10.1016/j.jneuroim.201 Neurology, Charité University Medicine Berlin, Berlin, Germany. 10Institute for 0.10.021. Computational Biomedicine, Heidelberg University Hospital and Faculty of 21. Saez-Rodriguez J, Alexopoulos LG, Epperlein J, Samaga R, Lauffenburger DA, Medicine, Heidelberg University, Bioquant, Heidelberg, Germany. Klamt S, et al. Discrete logic modelling as a means to link protein signalling networks with functional analysis of mammalian signal transduction. Mol Received: 2 May 2019 Accepted: 14 June 2021 Syst Biol. 2009;5(1):331. https://doi.org/10.1038/msb.2009.87. 22. Bernardo-Faura M, Massen S, Falk CS, Brady NR, Eils R. Data-derived modeling characterizes plasticity of MAPK signaling in melanoma. Plos Comput Biol. 2014;10(9):e1003795. https://doi.org/10.1371/journal.pcbi.1003795. References 23. Eduati F, Doldàn-Martelli V, Klinger B, Cokelaer T, Sieber A, Kogera F, et al. 1. Owens J. Determining druggability. Nat Rev Drug Discov. 2007;6(3):187. Drug resistance mechanisms in colorectal cancer dissected with cell type- https://doi.org/10.1038/nrd2275. specific dynamic logic models. Cancer Res. 2017;77(12):3364–75. https://doi. 2. Jia J, Zhu F, Ma X, Cao Z, Cao ZW, Li Y, et al. Mechanisms of drug org/10.1158/0008-5472.CAN-17-0078. combinations: interaction and network perspectives. Nat Rev Drug Discov. 24. Silverbush D, Grosskurth S, Wang D, Powell F, Gottgens B, Dry J, et al. Cell- 2009;8(2):111–28. https://doi.org/10.1038/nrd2683. specific computational modeling of the PIM pathway in acute myeloid 3. Cully M. Combinations with checkpoint inhibitors at wavefront of cancer leukemia. Cancer Res. 2017;77(4):827–38. https://doi.org/10.1158/0008-5472. immunotherapy. Nat Rev Drug Discov. 2015;14(6):374–5. https://doi.org/10.1 CAN-16-1578. 038/nrd4648. 25. Flobak Å, Baudot A, Remy E, Thommesen L, Thieffry D, Kuiper M, et al. 4. Klinger B, Sieber A, Fritsche-Guenther R, Witzel F, Berry L, Schumacher D, Discovery of drug synergies in gastric cancer cells predicted by logical et al. Network quantification of EGFR signaling unveils potential for targeted modeling. Plos Comput Biol. 2015;11(8):e1004426. https://doi.org/10.1371/ combination therapy. Mol Syst Biol. 2013;9(1):673. https://doi.org/10.1038/ journal.pcbi.1004426. msb.2013.29. 26. Kotelnikova E, Kiani NA, Messinis D, Pertsovskaya I, Pliaka V, Bernardo-Faura 5. Bozic I, Reiter JG, Allen B, Antal T, Chatterjee K, Shah P, et al. Evolutionary M, et al. MAPK pathway and B cells overactivation in multiple sclerosis dynamics of cancer in response to targeted combination therapy. Elife. revealed by phosphoproteomics and genomic analysis. Proc Natl Acad Sci 2013;2:e00747. https://doi.org/10.7554/eLife.00747. U S A. 2019;116(19):9671–6. https://doi.org/10.1073/pnas.1818347116. 6. Villoslada P, Steinman L, Baranzini SE. Systems biology and its application to 27. Bellmann-Strobl J, Paul F, Wuerfel J, Dörr J, Infante-Duarte C, Heidrich E, the understanding of neurological diseases. Ann Neurol. 2009;65(2):124–39. et al. Epigallocatechin gallate in relapsing-remitting multiple sclerosis. https://doi.org/10.1002/ana.21634. Neurol – Neuroimmunol Neuroinflamm. 2021;8(3) Available from: https:// 7. Saez-Rodriguez J, Blüthgen N. Personalized signaling models for doi.org/10.1212/NXI.0000000000000981. personalized treatments. Mol Syst Biol. 2020;16(1):e9042. https://doi.org/10.1 28. Saez-Rodriguez J, Simeoni L, Lindquist JA, Hemenway R, Bommhardt U, Arndt 5252/msb.20199042. B, et al. A logical model provides insights into T cell receptor signaling. Plos 8. Kholodenko BN, Hancock JF, Kolch W. Signalling ballet in space and time. Computat Biol. 2007;3:e163. https://doi.org/10.1371/journal.pcbi.0030163. Nat Rev Mol Cell Biol. 2010;11(6):414–26. https://doi.org/10.1038/nrm2901. 29. Kanehisa M. Toward understanding the origin and evolution of cellular 9. Kolch W, Halasz M, Granovskaya M, Kholodenko BN. The dynamic control of organisms. Protein Sci. 2019;28(11):1947–51. https://doi.org/10.1002/pro.3715. signal transduction networks in cancer cells. Nat Rev Cancer. 2015;15(9): 30. Gough NR, Adler EM, Foley JF. Cell signaling: Details, details, details. Sci 515–27. https://doi.org/10.1038/nrc3983. STKE. 2007;2007:eg9. https://doi.org/10.1126/stke.4072007eg9. 10. Davis MM, Tato CM, Furman D. Systems immunology: just getting started. 31. Martens M, Ammar A, Riutta A, Waagmeester A, Slenter DN, Hanspers K, Nat Immunol. 2017;18(7):725–32. https://doi.org/10.1038/ni.3768. et al. WikiPathways: connecting communities. Nucleic Acids Res. 2021; 11. Korkut A, Wang W, Demir E, Aksoy BA, Jing X, Molinelli EJ, et al. Perturbation 49(D1):D613–21. https://doi.org/10.1093/nar/gkaa1024. biology nominates upstream-downstream drug combinations in RAF 32. Davies M, Nowotka M, Papadatos G, Dedman N, Gaulton A, Atkinson F, et al. inhibitor resistant melanoma cells. Elife. 2015;4. https://doi.org/10.7554/eLife. ChEMBL web services: streamlining access to drug discovery data and 04640. utilities. Nucleic Acids Res. 2015;43(W1):W612–20. https://doi.org/10.1093/na 12. Palacios R, Goni J, Martinez-Forero I, Iranzo J, Sepulcre J, Melero I, et al. A r/gkv352. network analysis of the human T-cell activation gene network identifies 33. UniProt Consortium. UniProt: a worldwide hub of protein knowledge. Nucleic JAGGED1 as a therapeutic target for autoimmune diseases. Plos One. 2007; Acids Res. 2019;47(D1):D506–15. https://doi.org/10.1093/nar/gky1049. 2(11):e1222. https://doi.org/10.1371/journal.pone.0001222. 34. Meldal BHM, Forner-Martinez O, Costanzo MC, Dana J, Demeter J, 13. Lee MJ, Ye AS, Gardino AK, Heijink AM, Sorger PK, MacBeath G, et al. Dumousseau M, et al. The complex portal--an encyclopaedia of Sequential application of anticancer drugs enhances cell death by rewiring macromolecular complexes. Nucleic Acids Res. 2015;43(Database issue): apoptotic signaling networks. Cell. 2012;149(4):780–94. https://doi.org/10.101 D479–84. https://doi.org/10.1093/nar/gku975. 6/j.cell.2012.03.031. 35. Martinez-Pasamar S, Abad E, Moreno B, Velez de Mendizabal N, Martinez- 14. Ransohoff RM, Hafler DA, Lucchinetti CF. Multiple sclerosis-a quiet Forero I, Garcia-Ojalvo J, et al. Dynamic cross-regulation of antigen-specific revolution. Nat Rev Neurol. 2015;11(3):134–42. https://doi.org/10.1038/ effector and regulatory T cell subpopulations and microglia in brain nrneurol.2015.14. autoimmunity. BMC Syst Biol. 2013;7:34. 15. International Multiple Sclerosis Genetics Consortium, Wellcome Trust Case 36. Poussin C, Mathis C, Alexopoulos LG, Messinis DE, Dulize RHJ, Belcastro V, Control Consortium 2, Sawcer S, Hellenthal G, Pirinen M, CCA S, et al. et al. The species translation challenge-a systems biology perspective on Genetic risk and a primary role for cell-mediated immune mechanisms in human and rat bronchial epithelial cells. Sci Data. 2014;1(1):140009. https:// multiple sclerosis. Nature. 2011;476(7359):214–9. doi.org/10.1038/sdata.2014.9. 16. Bayat Mokhtari R, Homayouni TS, Baluch N, Morgatskaya E, Kumar S, Das B, 37. Terfve C, Cokelaer T, Henriques D, MacNamara A, Goncalves E, Morris MK, et al. Combination therapy in combating cancer. Oncotarget. 2017;8(23): Iersel M, Lauffenburger DA, Saez-Rodriguez J CellNOptR: a flexible toolkit to 38022–43. https://doi.org/10.18632/oncotarget.16723. train protein signaling networks to data using multiple logic formalisms. 17. Kotelnikova E, Bernardo-Faura M, Silberberg G, Kiani NA, Messinis D, Melas BMC Syst Biol. 2012;6:133, 1, doi: https://doi.org/10.1186/1752-0509-6-133. IN, et al. Signaling networks in MS: a systems-based approach to developing 38. Iorio F, Bernardo-Faura M, Gobbi A, Cokelaer T, Jurman G, Saez-Rodriguez J. new pharmacological therapies. Mult Scler. 2015;21(2):138–46. https://doi. Efficient randomization of biological networks while preserving functional org/10.1177/1352458514543339. characterization of individual nodes. BMC Bioinformatics. 2016;17(1):542. 18. Kieseier BC, Stüve O. Multiple sclerosis: combination therapy in MS--still a https://doi.org/10.1186/s12859-016-1402-1. valid strategy. Nat Rev Neurol. 2011;7(12):659–60. https://doi.org/10.1038/ 39. White BJ, Tarabishy S, Venna VR, Manwani B, Benashski S, McCullough LD, et al. nrneurol.2011.164. Protection from cerebral ischemia by inhibition of TGFβ-activated kinase. Exp 19. Conway D, Cohen JA. Combination therapy in multiple sclerosis. Lancet Neurol. 2012;237(1):238–45. https://doi.org/10.1016/j.expneurol.2012.05.019. Neurol. 2010;9(3):299–308. https://doi.org/10.1016/S1474-4422(10)70007-7. 40. Chiba K, Kataoka H, Seki N, Shimano K, Koyama M, Fukunari A, et al. Fingolimod (FTY720), sphingosine 1-phosphate receptor modulator, shows Bernardo-Faura et al. Genome Medicine (2021) 13:117 Page 16 of 16 superior efficacy as compared with interferon-β in mouse experimental 59. Tarcic G, Boguslavsky SK, Wakim J, Kiuchi T, Liu A, Reinitz F, et al. An autoimmune encephalomyelitis. Int Immunopharmacol. 2011;11(3):366–72. unbiased screen identifies DEP-1 tumor suppressor as a phosphatase https://doi.org/10.1016/j.intimp.2010.10.005. controlling EGFR endocytosis. Curr Biol. 2009;19(21):1788–98. https://doi. 41. Day GS, Rae-Grant A, Armstrong MJ, Pringsheim T, Cofield SS, Marrie RA. org/10.1016/j.cub.2009.09.048. Identifying priority outcomes that influence selection of disease-modifying 60. Pertsovskaya I, Abad E, Domedel-Puig N, Garcia-Ojalvo J, Villoslada P. therapies in MS. Neurol Clin Pract. 2018;8(3):179–85. https://doi.org/10.1212/ Transient oscillatory dynamics of interferon beta signaling in macrophages. CPJ.0000000000000449. BMC Syst Biol. 2013;7(1):59. https://doi.org/10.1186/1752-0509-7-59. 42. Kappos L, Wolinsky JS, Giovannoni G, Arnold DL, Wang Q, Bernasconi C, et al. 61. Schneider WM, Chevillotte MD, Rice CM. Interferon-stimulated genes: a Contribution of Relapse-Independent Progression vs Relapse-Associated complex web of host defenses. Annu Rev Immunol. 2014;32(1):513–45. Worsening to Overall Confirmed Disability Accumulation in Typical Relapsing https://doi.org/10.1146/annurev-immunol-032713-120231. Multiple Sclerosis in a Pooled Analysis of 2 Randomized Clinical Trials. JAMA 62. Takuwa Y, Okamoto Y, Yoshioka K, Takuwa N. Sphingosine-1-phosphate Neurol. 2020. https://doi.org/10.1001/jamaneurol.2020.1568. signaling in physiology and diseases. Biofactors. 2012;38(5):329–37. https:// 43. Canto E, Isobe N, Didonna A, MS-EPIC Study Group, Hauser SL, Oksenberg JR. doi.org/10.1002/biof.1030. Aberrant STAT phosphorylation signaling in peripheral blood mononuclear 63. Arish M, Alaidarous M, Ali R, Akhter Y, Rub A. Implication of sphingosine-1- cells from multiple sclerosis patients. J Neuroinflammation. 2018;15(1):72. phosphate signaling in diseases: molecular mechanism and therapeutic 44. Moreno B, Hevia H, Santamaria M, Sepulcre J, Muñoz J, García-Trevijano ER, strategies. J Recept Signal Transduct Res. 2017;37(5):437–46. https://doi. et al. Methylthioadenosine reverses brain autoimmune disease. Ann Neurol. org/10.1080/10799893.2017.1358282. 2006;60(3):323–34. https://doi.org/10.1002/ana.20895. 64. Rose DM, Alon R, Ginsberg MH. Integrin modulation and signaling in 45. Chen D, Ireland SJ, Remington G, Alvarez E, Racke MK, Greenberg B, et al. leukocyte adhesion and migration. Immunol Rev. 2007;218(1):126–34. CD40-mediated NF-κB activation in B cells is increased in multiple sclerosis https://doi.org/10.1111/j.1600-065X.2007.00536.x. and modulated by therapeutics. J Immunol. 2016;197(11):4257–65. https:// 65. Abram CL, Lowell CA. The ins and outs of leukocyte integrin signaling. doi.org/10.4049/jimmunol.1600782. Annu Rev Immunol. 2009;27(1):339–62. https://doi.org/10.1146/annurev. 46. Housley WJ, Fernandez SD, Vera K, Murikinati SR, Grutzendler J, Cuerdon N, immunol.021908.132554. et al. Genetic variants associated with autoimmunity drive NFκB signaling 66. Luo B-H, Carman CV, Springer TA. Structural basis of integrin regulation and and responses to inflammatory stimuli. Sci Transl Med. 2015;7(291):291ra93. signaling. Annu Rev Immunol. 2007;25(1):619–47. https://doi.org/10.1146/a 47. Hussman JP, Beecham AH, Schmidt M, Martin ER, McCauley JL, Vance JM, nnurev.immunol.25.022106.141618. et al. GWAS analysis implicates NF-κB-mediated induction of inflammatory T 67. Totzke J, Gurbani D, Raphemot R, Hughes PF, Bodoor K, Carlson DA, et al. Takinib, a cells in multiple sclerosis. Genes Immun. 2016;17(5):305–12. https://doi.org/1 selective TAK1 inhibitor, broadens the therapeutic efficacy of TNF-α inhibition for 0.1038/gene.2016.23. cancer and autoimmune disease. Cell Chem Biol. 2017;24(8):1029–39.e7. 48. Johnson BA, Wang J, Taylor EM, Caillier SJ, Herbert J, Khan OA, et al. 68. International Multiple Sclerosis Genetics Consortium. Multiple sclerosis Multiple sclerosis susceptibility alleles in African Americans. Genes Immun. genomic map implicates peripheral immune cells and microglia in 2010;11(4):343–50. https://doi.org/10.1038/gene.2009.81. susceptibility. Science. 2019;365(6460). https://doi.org/10.1126/science.aav7188. 49. Hassan NJ, Simmonds SJ, Clarkson NG, Hanrahan S, Puklavec MJ, Bomb M, 69. Bendall SC, Nolan GP, Roederer M, Chattopadhyay PK. A deep profiler’s et al. CD6 regulates T-cell responses through activation-dependent guide to cytometry. Trends Immunol. 2012;33(7):323–32. https://doi.org/10.1 recruitment of the positive regulator SLP-76. Mol Cell Biol. 2006;26(17):6727– 016/j.it.2012.02.010. 38. https://doi.org/10.1128/MCB.00688-06. 70. Terfve CDA, Wilkes EH, Casado P, Cutillas PR, Saez-Rodriguez J. Large-scale 50. International Multiple Sclerosis Genetics Consortium (IMSGC), Beecham AH, models of signal propagation in human cells derived from discovery Patsopoulos NA, Xifara DK, Davis MF, Kemppinen A, et al. Analysis of phosphoproteomic data. Nat Commun. 2015;6(1):8033. https://doi.org/10.1 immune-related loci identifies 48 new susceptibility variants for multiple 038/ncomms9033. sclerosis. Nat Genet. 2013;45(11):1353–60. 71. Bulusu KC, Guha R, Mason DJ, Lewis RPI, Muratov E, Kalantar Motamedi Y, et al. 51. Martinez-Forero I, Garcia-Munoz R, Martinez-Pasamar S, Inoges S, Lopez-Diaz Modelling of compound combination effects and applications to efficacy and de Cerio A, Palacios R, et al. IL-10 suppressor activity and ex vivo Tr1 cell toxicity: state-of-the-art, challenges and perspectives. Drug Discov Today. 2016 function are impaired in multiple sclerosis. Eur J Immunol. 2008;38(2):576– Feb;21(2):225–38. https://doi.org/10.1016/j.drudis.2015.09.003. 86. https://doi.org/10.1002/eji.200737271. 72. Halasz M, Kholodenko BN, Kolch W, Santra T. Integrating network 52. Beyer T, Busse M, Hristov K, Gurbiel S, Smida M, Haus U-U, et al. Integrating reconstruction with mechanistic modeling to predict cancer therapies. Sci signals from the T-cell receptor and the interleukin-2 receptor. Plos Comput Signal. 2016;9(455):ra114. Biol. 2011;7(8):e1002121. https://doi.org/10.1371/journal.pcbi.1002121. 73. Bernardo-Faura M, Rinas M, Wirbel J, Pertsovskaya I, Pliaka V, Messinis DE, et al. CombiMS source code Github; 2016. Available from: https://github. 53. Slavin-Chiorini DC, Catalfamo M, Kudo-Saito C, Hodge JW, Schlom J, com/saezlab/combiMS Sabzevari H. Amplification of the lytic potential of effector/memory CD8+ cells by vector-based enhancement of ICAM-1 (CD54) in target cells: implications for intratumoral vaccine therapy. Cancer Gene Ther. 2004; Publisher’s Note 11(10):665–80. https://doi.org/10.1038/sj.cgt.7700741. Springer Nature remains neutral with regard to jurisdictional claims in 54. Swaminathan G, Tsygankov AY. The Cbl family proteins: ring leaders in published maps and institutional affiliations. regulation of cell signaling. J Cell Physiol. 2006;209(1):21–43. https://doi. org/10.1002/jcp.20694. 55. Sanna S, Pitzalis M, Zoledziewska M, Zara I, Sidore C, Murru R, et al. Variants within the immunoregulatory CBLB gene are associated with multiple sclerosis. Nat Genet. 2010;42(6):495–7. https://doi.org/10.1038/ng.584. 56. Stürner KH, Borgmeyer U, Schulze C, Pless O, Martin R. A multiple sclerosis– associated variant of CBLB links genetic risk with type I IFN function. J Immunol. 2014;193(9):4439–47. https://doi.org/10.4049/jimmunol.1303077. 57. Sellebjerg F, Krakauer M, Khademi M, Olsson T, Sørensen PS. FOXP3, CBLB and ITCH gene expression and cytotoxic T lymphocyte antigen 4 expression on CD4(+) CD25(high) T cells in multiple sclerosis. Clin Exp Immunol. 2012; 170(2):149–55. https://doi.org/10.1111/j.1365-2249.2012.04654.x. 58. Galisteo ML, Dikic I, Batzer AG, Langdon WY, Schlessinger J. Tyrosine phosphorylation of the c-cbl proto-oncogene protein product and association with epidermal growth factor (EGF) receptor upon EGF stimulation. J Biol Chem. 1995;270(35):20242–5. https://doi.org/10.1074/jbc.2 70.35.20242.