Large-scale models of signal propagation in human cells derived from discovery phosphoproteomic data

Mass spectrometry is widely used to probe the proteome and its modifications in an untargeted manner, with unrivalled coverage. Applied to phosphoproteomics, it has tremendous potential to interrogate phospho-signalling and its therapeutic implications. However, this task is complicated by issues of undersampling of the phosphoproteome and challenges stemming from its high-content but low-sample-throughput nature. Hence, methods using such data to reconstruct signalling networks have been limited to restricted data sets and insights (for example, groups of kinases likely to be active in a sample). We propose a new method to handle high-content discovery phosphoproteomics data on perturbation by putting it in the context of kinase/phosphatase-substrate knowledge, from which we derive and train logic models. We show, on a data set obtained through perturbations of cancer cells with small-molecule inhibitors, that this method can study the targets and effects of kinase inhibitors, and reconcile insights obtained from multiple data sets, a common issue with these data.

the perturbation (red contour nodes) from the drug targets, until steady state. All edges are positive because we model perturbation flow, as evidenced by increases/decreases in phosphorylation, and not these signed changes themselves, and there are only two types of gates, OR/AND (i.e. the target is set to 1 if any/all of the input nodes is/are equal to 1). F. Scoring and generation best models. Models are scored based on the scores S ij for nodes predicted as perturbed or not perturbed, under the set of conditions analysed: S model =Σ nodes predicted P S ij -Σ DP nodes predicted C S i,j (i.e. true positive predictions (negative S ij ) + false positive predictions (positive S ij ) -false negative predictions (negative S ij )). The frequency of each edge in a subset of the best scoring models (represented here as the weight of edges in the background structure) is computed (numbers next to the example structures, e.g. for the intermediate node represented (red box), one of the two best models selected has the 2 leftmost edges combined with an 'AND' gate, and one has the leftmost gate only; for the sink example (yellow box) both best models have the left edge; for the integrator (blue box) both models have the "all edges'" option). G. Weights correction. Edges are "copied" in proportion of their frequency in the best models (with a maximum number called "cap"), increasing their sampling probability, i.e. information from the boxes at step d is combined with information from the boxes at step f. For intermediates (red boxes), the 'AND' sampling probability is the average of its latest value and its frequency in the latest best models. For the integrator example (blue box), if cap=2 the new possibilities are: [edge1 x 2, edge2, (edge1 + edge2 ) x 2, no edge], with respective probabilities [ 2/6 , 1/6 , 2/6 , 1/6 ] (edge1: 1 bin from step d + (50%*2=1) bin from step f; edge2: 1 bin from step d; AND with both edges: 1 bin from step d + (50%*2=1) bin from step f; no edge: 1 bin from step d). For the sink example (yellow box), the new possible inputs are: [edge1 x 3, edge2] with probabilities [0.75, 0.25] (edge1: 1 bin from step d + (100%*2=2) bins from step f; edge2: 1 bin from step d). For the intermediate example above, edge1and edge2 will be sampled with respective probabilities 3/5 and 2/5 (edge1: 1 bin from step d + (100%*2=2) bins from step f; edge2: 1 bin from step d + (50%*2=2) bins from step f), and the AND will be sampled with probability (0.5+0.5)/2 (average of the sampling weight in box d and the frequency in best models from box f). The numbers above each input represent the new sampling probabilities at the second generation, matching to step d (first generation). Figure 3. Proof of principle on reduced scope data. A,B. Sites with a perturbation pattern consistent with the canonical PI3K-AKT1-MTOR hierarchy are selected. C. The network built from these sites and drug targets is used for three independent optimisations. D. The three optimisations converged to similar scores (top panel) and stable averaged edge frequencies (lower panel). E. Averaged frequencies are used to extract a network of most likely paths from drug targets (green diamonds) to data (red hexagons) with various levels of tolerance around the highest frequency edge (10% represented). Panels represent the data for a given site under each of our three pairs of kinase inhibitors; subdivided panels= multiple peptides mapping to the site. F. Example data and GMMs used to compute the scores (blue=control distribution, yellow=perturbed distribution). Note: all node names are Uniprot identifiers (without the _HUMAN suffix).  E. cap F. direct interactions parameters tested. G1=first generation of optimisation. Scores are shown without any size penalty because we found model size to be almost identical under any combination of parameters, due to the way that the sampling works. Unless stated otherwise, all optimisations are performed with settings direct interactions=G1, cap=5, n=5000, tolerance=30%, sizeP=1. A. The relative tolerance is the window of scores around the best score of each generation that is used to select the best models (in % of the best score). This parameter seems to have little influence on the scores between 15 and 50%, performing worse at the best model level at 70%. B. n is the number of models sampled at each generation. Unsurprisingly, higher n leads to better best models. In terms of average population scores however, most of the curves stabilise around the same area, with the lower numbers of models being more prone to outliers on either side (i.e. better or worse scores). In practice we found 500 to be a good compromise between optimisation time and stability for most models tested. C. The absolute tolerance is the same as the relative tolerance except that the window is fixed throughout generations. There is here a clear gradient from better (low tolerances) to worse (high tolerances). Based on this much higher sensitivity and difficulty to set a score range a priori, we would suggest to use relative tolerances. D. The size penalty is a factor that, multiplied by the fraction of edges sampled in each model, can be added to the model scores to account for model size in the optimisation process. At the population average level, the higher the size penalty, the bigger the variability between independent optimisations. In terms of best model scores, the highest size penalty shows poor performance. Size penalties have close to no influence on the size of models (due to the nature of the sampling process, most nodes will pick an incoming edge, so it is not so much the number of edges as their identity that varies from one model to another). However, higher size penalties allow poorer models with slight differences in size to sometimes reach better scores, resulting in a bigger variability between optimisations. We therefore excluded the size term in all optimisations. E. The cap is the maximum number of times that an edge is "copied through" in the background network, at the end of a generation (see Supplementary Fig. 2). This parameter determines the speed of convergence by influencing how much of an advantage edges in best models have at the next generation. The higher the cap, the better the average population score. The best model scores however seem to be largely unaffected. F. The "direct interactions" setting refers to the possibility to force include the direct (non-predicted) interactions between sites perturbed and drug targets at all generations, the first generation or not at all (in order to simplify the training process by forcing obvious interactions in or giving them an advantage in starting weights). Clearly, including the direct interactions at all generations produces better scores and allows the optimisation to focus on capturing other sites.  inhibition (A including, and B excluding predicted edges). A. All paths are represented at the protein level, for clarity (apart from the step between perturbed sites and all their known potential kinases). In consequence, every edge can represent multiple underlying edges, one for each site on the target K/P that is phosphorylated/dephosphorylated by the source K/P. Clearly, without the optimisation steps in PHONEMeS, this knowledge is of very limited use to interpret the data. B. Excluding predicted interactions, the network becomes easier to manipulate (261 edges and 163 nodes) but can only capture 15 of the 32 sites perturbed under MTOR inhibition. To capture the remaining 17 sites in our data, in the absence of a method like PHONEMeS that can pick out the kinases likely to lie on perturbation paths and connect predicted sites to these kinases, the only option available is to use the full (hardly visualisable) background in A. . Most sites are connected to the designated targets in all cases but the paths are shorter with the correct target than with the incorrect ones.

Supplementary
The random networks show a major difference in terms of number of sites explained directly by the drug target or its first neighbours (14 sites for the real network vs. 4 sites for this random network). Figure 8. Modelling the effect of MTOR inhibition in a publicly available data set 13 . A. 12 sites perturbed under one of the two MTOR inhibitors in the validation data were connected to MTOR in six independent optimisations, using exclusively experimentally demonstrated interactions. Results of these were combined into average frequencies in final populations of models, which were used to extract a best consensus model with 5% tolerance around the highest frequency edges. This model recapitulates effects on known substrates of MTOR, and shows a good agreement with the data obtained with a different mass spectrometer and reagents in Fig. 4. This network is less constrained and shows a lower agreement between the two drugs when looking at sites further downstream from MTOR. This could be due to the lower coverage, and possible lower quality of this data set (due partially to the type of spectrometer used, and to the data processing approach, see Wilkes et al. 22 ). B. Sites that were found to be perturbed but did not have an experimentally demonstrated kinase were connected to the optimised network based on predictions from networKIN 17 . Note that sites which could be connected to a kinase that is a first neighbour of mTOR are more likely to be consistently identified under both drug treatments, further supporting our network-based prioritisation of potential kinases. C. Predictions were prioritised based on the sign of the perturbations observed, and functional annotations of proteins were extracted from Uniprot.

Supplementary
We can see that many of the functional themes and network connections that were identified with the main data in Fig. 4 are also found here, even though the two sets of experiments involved different drugs, spectrometric equipments and as a consequence vastly different data sets and perturbed peptides. Additionally, this analysis shows that even with data with lower coverage and less advanced technology it is possible to obtain a network that is informative both in terms of the functional effect of a drug treatment and in terms of the prioritisation of kinase-substrate relationships downstream of this drug treatment. Figure 9. Inhibition of mTORC1/2 or Akt decreases phosphorylation of Cyclin-L1. A. Extracted ion chromatograms for the ion representing two phosphorylation sites on Cyclin-L1 (CCNL1) under each of the inhibitor treatments (blue lines, first isotope; red lines, second isotope; green lines, third isotope). B. Box plots summarising the raw intensity data shown in (B). Dotted line represents the median intensity under treatment with DMSO vehicle control. This phosphorylation of this CCNL1 peptide was decreased 2 to 5 fold following Akt and mTORC1/2 inhibition, as well as PI3K inhibition, but not P70S6K, which is consistent with the network prediction in Fig. 4. size (number of edges) and scores along 50 generations of optimisation when both drugs are assumed to target AKT1 and AKT2 (red, green and yellow curves, AKTib) or when AKTi2 is assumed to also target AKT3 (blue and purple curves). Bottom: network resulting from the combined output of the 3 rounds of optimisations AKTib, visualised with 0% tolerance around the highest frequency input edges. Edges widths represent the average frequency at which the edge is present in the population of networks at the last generation of the optimisation.

Supplementary
Supplementary Figure 11. Modelling the effect of the MEK inhibitors GSK-1120212 (MEKi1) and U0126 (MEKi2). Network resulting from the combined output of three rounds of optimisations with both drugs assumed to have the same targets MP2K1 and MP2K2, visualised with 20% tolerance around the highest frequency input edges. Edges widths represent the average frequency at which the edge is present in the population of networks at the last generation of the optimisation.

Supplementary Figure 12. Modelling the individual effect of the Erk inhibitors I (ERKi2) and II (ERKi2).
Results of three rounds of optimisation with different potential targets for ERKi1 and ERKi2. Top: networks resulting from the combined output of the three rounds of optimisations ERKi1b (ERKi1 with targets MK01 and MK03) and ERKi2b (ERKi2 with targets MK01 and MK03), visualised with 10% and 0% tolerance around the highest frequency input edges, respectively. Edges widths represent the average frequency at which the edge is present in the population of networks at the last generation of the optimisation. Bottom: population average (continuous) and best model (dashed) scores along 50 generations of optimisation: left, when ERKi1 is assumed to target MK01 and MK03 (red, light green and orange curves), MK01, MK03 and MK14 (light blue, turquoise and dark green curves) or MK01 only (dark blue, pink and purple curves); right, when ERKi2 is assumed to target MK03 only (red, orange and yellow curves), MK01 only (blue curves), MK01 and MK03 (green curves) or MK01, MK03 and MK14 (pink and purple curves).  Figure 16. Combined modelling of the effect of the P70S6K inhibitors PF-4708671 (P70SKi1) and S6K Inhibitor (P70S6Ki2). Results of three rounds of optimisation with different potential targets for P70S6Ki1 and P70S6Ki2. Top: network resulting from the combined output of three rounds of optimisation in which both drugs are assumed to have KS6B1 as their only target, visualised with 20% tolerance around the highest frequency input edges. Bottom: population average (continuous) and best model (dashed) scores along 50 generations of optimisation: left, for P70S6Ki1 and right, for P70S6Ki2, when they are both assumed to have KS6B1 as their only target (pink and blue curves) or when P70S6Ki1 is assumed to additionally target KS6B2. Network and data representation conventions are as in Supplementary Figs.10-12.

Supplementary Figure 13. Combined modelling of the effect of the Erk inhibitors I (ERKi2) and II (ERKi2
Supplementary Figure 17. Analysis of the mTOR inhibition data by Hsu et al. 13 . Published iTRAQ labelled data obtained after inhibition of mTORC1 by Rapamycin (B) and mTORC1 and mTORC2 by Torin1 (B) in HEK-293E cells (in insulin-stimulated cells) was analysed using PHONEMeS. The networks show the results of PHONEMeS analyses run separately for the two drugs, shown with 20% tolerance around the best incoming edges. The data shows the cluster assignment (pink/blue=more than 2.5MAD larger/smaller than the average log 2 (drug/insulin)), and the average log 2 (drug/insulin) for the peptide mapping to the indicated phosphorylation site (two replicates). Red nodes are sites that are perturbed by the drug treatment, grey nodes are nodes for which no information is available in the data, green nodes are nodes that are designated as drug targets. The edge widths reflect the frequency with which the edge is present in the population of models at the end of the optimisation.   A GMM is fitted to each peptide and those with two components are kept. Scores Si,j are computed for each data point (peptide i, drug j) that capture the log ratio of the probability of being in the control vs perturbed distributions, using parameters estimates from the GMM. If a measurement has an S ij below -0.5 (i.e. it is highly more likely to be in the perturbed state) and the FC vs control has an adjusted p-value below 0.05 (i.e. the effect of the treatment is reproducible across replicates), then we refer to it as a "double positive". If a measurement has an S ij above 0.5 and an adjusted p-value for FC vs control above 0.05, then we call it a "double negative" (see Supplementary Table 3). The two metrics (FC significance and score) agree in most cases, with 78.6% of the data being correctly categorised as either double negatives or positives (DP, DN). A small proportion of peptides lie just in between the control and perturbed distribution (U, i.e. undetermined), and a small proportion of significantly perturbed peptides are misclassified as control (SN, single negative, i.e. significant FC classified in the control distribution). A larger proportion is classified as single positives (SP), which is expected with MS data as a result of relatively high FCs with low reliability (e.g. when one of the data points, either the control or the drug treated point, lie at an extreme of the dynamic range of the machine).
S i,j < -0.5 -0.5 < S i,j < 0.  This table contains the kinases predicted by NetworKIN as potentially responsible for the phosphorylation of sites found to be perturbed under MTOR inhibition, and for which no known kinase was available in other databases. By connecting these sites to the predicted kinases in our optimized network, and combining with experimental evidence of perturbation, we were allowed to prioritise 1 (or 2) kinases that are likely to be responsible for the phosphorylation of these sites, at least in the conditions observed.

Supplementary Notes
Supplementary Note 1

Proof of principle on data with reduced scope
In order to test our method, we reasoned that it should be able to recover a canonical kinase hierarchy when exposed to a reduced dataset that is consistent with this hierarchy. Therefore, we built and trained models based on the data for 6 small-molecule inhibitors (in pairs targeting class I PI3-kinases, Akt and mTOR) and 19 sites whose behaviours are consistent with the canonical view that PI3K leads to the activation of Akt which is then responsible for the activation of mTOR ( Supplementary Fig. 3a-b, f), such that sites affected by inhibition of mTOR or Akt should also be affected by inhibition of PI3K etc. Twenty-three additional sites that lie on potential paths between the 19 perturbed sites and the drug targets were included in the problem at the automatic background network-building step ( Supplementary Fig. 3c). We also used this setting to study the influence of all free training parameters in order to determine the reasonable ranges to use (Supplementary Fig. 4). This approach does not impose any constraint on noise, function or identity of the sites selected. It does however allow for the presence of contradictory and uncertain information on paths between targets and data sites, which would unavoidably occur in a real setting.
Three independent trainings were combined into a single solution ( Supplementary Fig. 3d). The resulting network contained 39 nodes (all referred to by their Uniprot identifiers) and 36 edges with 0 % tolerance (12 and 6 % of the total nodes and edges in the background network, respectively), and 45 nodes and 49 edges (15 and 8 %) with 20 % tolerance ( Supplementary  Fig. 3e). This network robustly captures the expected hierarchy of inhibited kinases (even with noisy and inconsistent data) and the perturbation of many of their canonical substrates (Fig.  3e). Indeed, the solutions correctly recovers the PK3CA->AKT1->MTOR cascade and connects the majority or sites to these or downstream kinases, including sites that are not direct targets of any of the inhibited kinases. It also captures perturbations of "unexpected" sign (i.e. phosphorylation increased upon kinase inhibition), such as the increase in DPYL2-pT509, which is explained in the model by the decrease in phosphorylation of the inhibitory GSK3B-pS9. The model cannot fully capture some sites found perturbed under PI3K inhibition (SRPK1-pS51, HNRPC-pS253 and possibly CDK16-pS153 and SRRM1-pS769), which is an expected consequence of PI3K being primarily a lipid kinase and not a protein kinase (whereas our background network only covers protein substrates). Within these limitations, PHONEMeS has demonstrated its ability to recover a canonical kinase signaling hierarchy when the data comply with its implied signal flow. Therefore, using unfiltered phospho-MS data with PHONEMeS should allow us to study non-canonical aspects of kinase inhibitions and resulting downstream perturbations.

Single drugs analysis: efficacy and specificity in context
The data set analysed here contains drug treatments in pairs with the same nominal target(s), to increase confidence in the data by focusing on perturbations observed under both drugs in a pair (hence reducing the effects of noise and off target effects). However, there are two reasons for a drug pair to show a different set of perturbed sites: i) the two drugs have slightly different targets (or off target) effects, ii) a slightly different set of sites under the same targets are detected and/or reach significance. Secondary targets are common and often poorly characterised, especially when studying kinase inhibition in cells rather than in solution. However, when looking at shotgun MS data, option ii) is also a plausible explanation for differences between drugs. Using PHONEMeS, we can look at the efficacy and specificity of drugs in a pair in the presence of both of these sources of uncertainty. In order to do this, we separately optimised a network for each drug treatment and each pair of drug treatments, using different sets of potential targets based on the drugs' nominal targets and secondary targets mentioned by the manufacturer, and secondary targets found in chEMBL (see Supplementary Table 7).
Because additional targets can make the background network significantly more complex (see table Supplementary Table 7), leading to less stable and less determined solutions, we decided to be relatively conservative in adding secondary targets. In addition, in some cases adding a target to one drug in a pair will allow to explain sites that are different between two drugs in a pair, regardless of whether they result from actual differences in targets or just in detection of the perturbations. However, if there really is a difference in drug specificity, then adding the target to the individual drug optimisation should also lead to significant score improvements. Additionally, when the candidate target is already pulled in a network where it is not specified as a target (i.e. it seems necessary to explain the data) we can be more confident that the candidate kinase is a "real" target of the drug. Finally, when interpreting the results of this analysis, it is important to remember that drugs can show widely variable effects on the network both in terms of spread (how far they reach) and breadth (how many areas of the network they hit) for multiple reasons: i) if the kinase inhibited is highly specific to a limited number of substrates, the breath of the perturbation will accordingly be small; ii) we are looking at a single snapshot in time, hence when the signalling spreads slowly in the network, it is possible that the full reach of an inhibition has not yet been realised; iii) drugs can have partial inhibitory activity, which is difficult to capture in a Boolean setting; iv) we are looking exclusively at inhibition compared to the basal state, hence we will not see anything in the network if the kinase was not significantly activated in the basal state.
All optimisations were performed as described in the main text. The results of all of these analyses can be found in table Supplementary Table 7 and the associated resulting networks can be downloaded as Cytoscape workspaces on the PHONEMeS webpage http://www.cellnopt.org/phonemes. Some of these analyses are discussed in detail below.

Akt
Akt (or protein kinase B) is a serine/threonine kinase that regulates key cellular processes from cell growth to survival. Mammalian cells express three different isoforms of Akt (Akt1,2,3) with close specificity but different tissue-specific patterns of expression 7,11 . Akt1 activation by growth factors is relatively well known and relies on the membrane recruitment of Akt via its pleckstrin homology (PH) domain which interacts with phosphoinositides PI(3,4)P2 and PI (3,4,5)P3 produced by activated phosphatidylinositol 3-kinase (PI3K). This interaction enables phosphorylation of Akt in its activation loop at T308 by PDPK1. mTORC2 or DNA-PK phosphorylation of Akt at S473 results in a fully activated Akt 7,11 . Activated Akt translocates to the cytosol and nucleus, where it phosphorylates many substrates.
AKTi1 is a PH domain-dependent inhibitor of Akt1,2 (and Akt3 to a lesser extent), and AKTi2 is an allosteric non-ATP competitive inhibitor of Akt1,2 and 3. Both inhibitors are, according to manufacturer information, specific inhibitors of these targets, which their chEMBL records seem to confirm. However, AKTi1 treatment results in many more sites perturbed than the AKTi2 treatment. As we can see from table Supplementary Table 7 and Supplementary Fig.  10, the score (for both drugs) are much better when Akt1 and Akt2 are used as common targets of both drugs (final population average scores improved from 11.3±1.3 to 7.2±1.5 for AKTi1 and 5.7±0.3 to 3.6±0.9 for AKTi2), and the models are also much smaller when Akt3 is not included (final population average model size from 892/892/894 to 591/592/615 edges). Additionally, even when Akt3 is included in the AKTi2 optimisation, few sites (2 out of 71) are actually placed under Akt3 (not shown), indicating that Akt3 is probably not necessary to explain the data. This would suggest that either Akt3 is not significantly inhibited by the drugs, or it is simply not significantly active (or even expressed) in this cell line under these conditions. Finally, in the combined optimisations, we did not see sites perturbed under only one of the two drugs clustering under any specific kinase, indicating that the difference between the two drugs are more likely to result from a difference in detection rather than specificity (see Supplementary Fig. 10). The fact that, for a lot of sites perturbed only under AKTi1, the sign of the fold change under AKTi2 is consistent with the one under AKTi1 but just doesn't reach significance, supports this hypothesis (see Supplementary Fig. 10). Together, these results suggest that both AKT inhibitors are best represented by inhibition of Akt1 and Akt2.

MEK
MAPK/ERK kinase (MEK) 1 and 2 are two closely related dual specificity protein kinases, part of the mitogen-activated protein kinase (MAPK) cascades which are involved in regulation of cell proliferation, survival and differentiation. MEK 1/2 are activated by the Raf serine/threonine kinases (c-Raf-1, A-Raf, B-Raf), which are themselves activated primarily by Ras small GTPases (although additional signalling activities through PAK, Src, PP2A, 14-3-3 proteins etc have been shown to play a role in Raf activation). The Ras proteins (H-, K-and N-Ras) are activated by extracellular stimuli such as EGF through the epidermal growth factor receptor (EGFR). MEK1 and 2 phosphorylate serine/threonine and tyrosine residues of their only known substrates ERK1/2 (extracellular signal-regulated kinase 1/2). Activated ERKs regulate a large number of (mainly nuclear) proteins. Mutations and/or overexpression at multiple levels of this cascade are frequently implicated in a variety of cancers 21 .
MEKi1 is a selective allosteric inhibitor of Mek1,2 kinase activity, and MEKi2 is a specific noncompetitive (with respect to ATP and Erk) inhibitor of Mek1,2. As we can see from Supplementary Table 7, these drug treatments result in the lowest number of perturbed sites of all drugs tested. Given the narrow substrate specificity of the MEK kinases and the high specificity of their inhibitors 21 , it is not overly surprising that the number of sites perturbed is relatively small, although we would have perhaps expected to see more of the downstream effects via substrates of ERK. It is possible that the time frame used in this experiment did not allow for this to appear quite strongly enough to be consistently detected. Nonetheless, the reconstructed paths found here accurately reflect the MEK-ERK-JNK expected signalling flow and show defined paths to plausible effectors of the pro-survival activity of this pathway. Both drugs have the expected result of a decrease in phosphorylation of MK01-Y187 and MK03-Y204, which are activating sites of Erk1 and Erk2, respectively (see Supplementary Fig. 11.). Paths of interest include the decrease in phosphorylation of B2L11-S77, which is an inducer of apoptosis, under MK10-T221 (which is an activating site of MK10/JNK3, and the decrease of B2L11-S77 is therefore consistent with this path). This site is not recorded as functional in Uniprot but interestingly S69 is reported to be phosphorylated by Erk1,2, leading to proteasomal degradation of B2L11. The decrease in phosphorylation of TR10A.S466 (reached via MK03 and KC1E) is also of interest since this protein is the receptor for the cytotoxic ligand TRAIL, and therefore also plays a role in induction of apoptosis.

ERK
Erk1/2 are terminal serine/threonine kinases of the MAPK cascades activated by growth factor-stimulated cell surface receptors, such as receptor tyrosine kinases (RTKs) of EGF and insulin 2,21 . Upon activation by MEK1/2 (see above), ERKs can translocate to the nucleus where they regulate a variety of transcription factors, leading to changes in cell proliferation, survival and motility. Erk1/2 signalling can also regulate cytoplasmic substrates like TSC2, RSK, members of the BCL2 family of apoptotic regulators, etc. Erk1/2 also exert negative feedback effects at multiple levels in the MAPK cascade 2, 18 .
ERKi1 is a non-ATP-competitive erk1/2 inhibitor that prevents interaction with substrates by binding erk docking domains. It preferentially binds erk2, and has little effect on erk1/2 phosphorylation by mek1/2 10 . ERKi2 is a potent ATP-competitive inhibitor of erk1/2 (MK01/MK03), which according to chEMBL also has activity towards MK14. As we can see in Supplementary Table 7, the two drugs only have four perturbed sites in common. Considering that the two drugs have different modes of action and partially overlapping known targets, this is perhaps not overly surprising but constitutes a sure warning sign. While performing the separate optimisations for ERKi1 and ERKi2, we noticed that the ERKi1 optimisation performed best when MK14 is added as an off-target of the drug (final population average score from 0.52±0.02 to 0.40±0.07, see Supplementary Table 7. Additionally, MK14 was pulled into the solution networks when it was not specified as a target, with direct substrates showing a decreased phosphorylation level (see Supplementary Fig. 12). On the other hand, the ERKi2 optimisations showed no improvement upon addition of MK14 as a target (final population average score from 0.69±0.22 to 0.67±0.27, and MK14 was not pulled into the network when it was not specified as a target. Instead, it is the optimisation with MK03 only that performs best (final population average score of 0.24±0.02, see Supplementary Fig. 12 and Supplementary Table 7. It would therefore seem that the identity of the two drugs was swapped somewhere in the process of generating or recording the data. As we can see on Supplementary Fig. 13, optimising with MK01, MK03 and MK14 as targets of ERKi1, and MK03 as a target of ERKi2 performed best in terms of scores (final population average score of 1.1±0.2 compared to 1.5±0.1, 1.7±0.4 and 1.7±0.1 for the other setups, see Supplementary  Table 7 but also allows placing of perturbed sites in a much more consistent way (compare the networks on the left and on the right on Supplementary Fig. 13: sites found to be perturbed only under ERKi2 are mostly placed under MK03, whereas many sites found to be perturbed exclusively under ERKi1 are placed under MK14). This setup manages to relatively successfully explain some of the large differences in effect of the two drugs that we have seen above.

PI3K
The PI3Ks are lipid kinases activated by growth factor receptor tyrosine kinases (RTK) and G-protein-coupled receptors (GPCRs). PI3Ks phosphorylate phosphatidylinositol-4,5-bisphosphate (PIP2) to produce phosphatidylinositol-3,4,5-trisphosphate (PIP3). Phosphatidylinositol is a membrane phospholipid that can be phosphorylated at the 3, 4 and 5 positions of the inositol ring, generating species with different interaction properties and signalling outcomes (see Bunney et al. 1 for a review). PIP3 recruits to the plasma membrane and activates several PH domain-containing proteins such as PDK1 and AKT, driving cell cycle progression and survival 5,8,19 . Class IA PI3Ks are heterodimeric lipid kinases containing a p110 catalytic subunit (PIK3CA, PIK3CB and PIK3CD) and a p85 regulatory subunit. PIK3CA mutations are the most common genetic alterations of this pathway in breast cancer (Miller et al., 2011), and indeed the cell line used here, MCF7, harbours a E542K mutation within the helical domain of p110α ( s e e http://cancer.sanger.ac.uk/cosmic/mutation/overview?id=760).
Such mutations seem to abrogate the inhibitory interaction between p85 and p110, hence the mutant p110α subunits are known to increase in vitro lipid kinase activity and, in some cases, to maintain PI3K-Akt signalling under growth-factor deprivation 5 . Cancer cells expressing those mutants are however similarly susceptible to PI3K inhibitors as those expressing wild type PI3K 3 , and studies in patients seem to indicate that such mutations confer increased sensitivity to PI3K/Akt/mTOR inhibitions, possibly indicating a dependency of these cancers on this pathway 14 .
PI3Ki1 is a potent ATP-competitive pan-PI3K inhibitor (equally potent against oncogenic mutants of p110α) 19 . On top of the nominal targets (PK3CA, PK3CB, PK3CG, PK3CD, chEMBL records AKT1 and MTOR as targets of this inhibitor. PI3Ki2 is an ATP-competitive inhibitor of PI3K, mTOR and DNA-PK. Small molecule inhibitors of p110 often also inhibit mTOR due to structural similarities. The chEMBL record for PI3Ki2 only lists interaction with AKT1. When optimising with PI3Ki1 alone, the setup with AKT1 and MTOR as off targets produced the best scores (final population average score of 0.88±0.04 vs 1.72±0.19, see Supplementary  Table 7) and the decrease of phosphorylation at many direct targets of these two kinases would confirm that they are affected by the drug (see Supplementary Fig. 14, left). Since these two kinases are the main downstream targets of activated PI3K, it is also possible that the improvement in the optimisation comes from the fact that declaring those as targets simplifies the search by automatically connecting their perturbed substrates. A similar situation is observed for PI3Ki2, where adding AKTi1 as an off target drastically improved the scores from final population average scores of 4.45±0.17 to 1.42±0.47, see table ST.7. MTOR was pulled into the network whether it was declared as a target or not, and adding it as a target also improved the scores to 2.10±0.29 (see Supplementary Fig. 14, right and Supplementary Table 7). In comparison, specifying AKT1 as a target did improve the scores (see above), but it seemed to do so mainly by explaining its direct substrates. This would indicate that both drugs have PI3K (PK3CD and PK3CA) and MTOR as main targets. The combined network scores were drastically improved by adding off targets MTOR and AKT1, especially for PI3Ki2 (final population average scores from 1.82±0.07 to 1.01±0.32 for PI3Ki1 and 4.16±0.05 to 1.38±0.27 for PI3Ki2, see Supplementary Table 7. The precise connectivity between MTOR, AKT, PI3K and KS6B1 is somewhat unclear because most configurations are equivalent as far as the model is concerned, but there is clearly an effect of both drugs on MTOR and AKT1 (either directly, which seems likely for MTOR, or via PI3K/MTOR inhibition, which seems more likely for AKT). Beyond the clear effect on MTOR and AKT, the networks suggest an important role for the kinase SRC.

P70S6K
The 70kDa ribosomal protein S6 kinases (S6K1 and S6K2, collectively P70S6K) are serine/threonine kinases of the AGC family (which also includes Akt, PKCs, P90RSK, etc). The catalytic domains of S6K1 and S6K2 are very similar but differences in the extreme C-and N-terminal regions are thought to lead to different localisation and target properties. The S6Ks are primarily activated through the PI3K pathway. The physiological target of S6Ks is the ribosomal S6 protein (rs6p, a component of the 40s ribosomal subunit), the phosphorylation of which plays a key role in modulating translational efficiency. S6Ks also directly phosphorylate initiation (e.g. eiF4B) and elongation (e.g. eEF2k) factors, as well as other proteins involved in proliferation and survival. Additionally, S6Ks seem to play an important role in mediating negative feedback loops in the PI3K signalling network 6 .
P70S6Ki1 is an ATP-competitive S6K1-specific inhibitor that was shown to prevent S6K1mediated phosphorylation of rs6p. Both the single drugs and combined drugs network (see Supplementary Fig. 15 and Supplementary Fig. 16) look cleaner when only KS6B1 (S6K1) is used as a unique common target of the two drugs, although the score curves are remarkably variable (final population average scores of 1.46±0.76 and 1.34±0.20 for the P70S6Ki1 optimisations, 2.35±0.45 and 2.70±1.05 for the combined optimisations, respectively with and without KS6B2, see Supplementary Table 7). This makes it hard to conclude anything on the scores differences, especially since the problem is worsened by a doubling of the size of the networks after adding S6K2 as a target (background network from 853 to 2292 edges for P70S6Ki1 optimisations, and 1133 to 2809 edges for the combined optimisations, see Supplementary Table 7).
Despite differences in sites perturbed, many paths found downstream of P70S6Ki1 are also found under P70S6Ki2 (compare Supplementary Fig. 15 top vs. bottom), and there does not seem to be a clear clustering of sites specific to one or the other drug under specific kinases in the combined network ( Supplementary Fig. 15), further indicating that the two drugs are likely to have the same target. Unfortunately, beyond the direct effect on RICTR and the apparent feedback on PDK1 (via an increase in its autophosphorylated activation site S241) on which the two drugs agree, the paths observed are fairly indirect and unstable. This seems to be a consequence of the fact that many of these sites are connected to predicted rather than known kinases (which are less specific links and hence typically less stable in the optimisation). It is therefore possible that, in addition to a far from perfect agreement between the data related to the two drugs, these networks also suffer from a relative lack of precise knowledge about substrates of S6K1 (see Fenton et al. 6 for a review).

Analysis of a published unrelated MS data set
In order to prove the general applicability of our method, we went on to analyse a published data set 13 . This data set consists of iTRAQ labelled phosphoproteomic samples after 1 hour treatment of HEK-293E cells with 100 nM rapamycin, 250 nM Torin1, or DMSO followed by 20 minutes stimulation with 150 nM insulin. The data (supplementary table S2 of Hsu et al. 13 ) consists of drug treatment vs insulin log2 ratios in biological duplicates, across 3883 peptides of which only the 1609 having data for both replicates were kept. This data was centred and standardised by replicate, and the replicate measurements were averaged for each peptide. The Uniprot identifiers of these peptides were obtained by mapping the RefSeq identifiers provided using the Uniprot ID mapping tool (http://www.uniprot.org/uploadlists/). This data is much less suited for quantitative modelling of phosphorylation networks than the data that we generated ourselves. The limited number of replicates available means that it is not possible for us to obtain statistically meaningful estimates of the effect of the treatments on each peptide. Additionally, because only two treatments are applied and the data is by nature in the form of ratios rather than absolute values for each treatment and a matched control, we cannot estimate a perturbed and a control state for each phosphorylation site with a Gaussian mixture model, as we would ideally do. However, the scores used for the optimisation step are derived from the Gaussian mixture model estimates (i.e. they capture how much more likely a peptide is to be in the perturbed rather than the control state under a certain condition). We therefore had to adapt the scoring procedure in order to assign to each peptide under each of the two treatments a score that is negative if the peptide is more likely to be perturbed, and which decreases as the amplitude of the difference from the control increases. In their analysis, Hsu et al. use peptides that are 2.5 median absolute deviation below the average of the condition to define the set of peptides that are affected by a drug (under the assumption that most peptides are not expected to react to a treatment when using high-content untargeted data). As a consequence, they are essentially selecting the peptides that show a large decrease in phosphorylation under torin or rapamycin treatment, compared to the insulin (DMSO) control.
We are however interested in peptides that show an increase in phosphorylation upon treatment because our method can (and aims to) capture more than direct targets of mTOR. Nonetheless, we wanted to stay close to the strategy used by Hsu et al. 13 , in order to be able to put our results in perspective with theirs. Therefore, we decided to use the following scoring scheme: the score of a peptide under a drug treatment (rapamycin or mTOR compared to insulin control) is -log 10 (p/1-p) where p is the probability of obtaining a value further from the mean of the distribution for that drug treatment. As required for our optimisation scheme, this value is negative and its absolute value increases as the effect observed moves further from the mean. However, whereas our score can be negative (if a data point is more likely to belong to the perturbed distribution) or positive (if a data point is more likely to belong to the control distribution), this score is always negative because the above ration is always smaller or equal to 1. This problem can be handled by the other element that is used in our optimisation scheme, which is the assignment of each data point to the control or perturbed cluster. In order to stay close to the Hsu et al. 13 analysis, we assigned the "P" or "C" status based on a 2.5 MAD distance from the average of the condition. The status of this assignment is further called single negative/single positive if the two replicate measurements were too far from each other, in order to exclude clearly unreliable peptides, in the absence of a proper measure of significance. This is done by filtering the peptides using a threshold on the sum of the differences between the replicates in a pair divided by their mean. A threshold of 10 (arbitrarily chosen to exclude the most unreliable peptides while preserving a reasonable amount of data) led to a total of 1042 peptides (i.e. 567 peptides were classified as single negatives/positives).
As mentioned above, mTOR functions as part of two different complexes, mTORC1 and mTORC2. mTORC2 signalling is insensitive to nutrients but responds to growth factors such as insulin via PI3K. mTORC2 directly activates Akt by S473 phosphorylation, as well as serum-and glucocorticoid-induced protein kinase 1 (SGK1, which controls ion transport and growth), and protein kinase C α (PKCα, which is involved in regulation of cell shape) (Laplante and Sabatini, 2012). mTORC1 integrates inputs from growth factors, stress, energy status, oxygen and amino acids. The heterodimer tuberous sclerosis 1 and 2 (TSC1 and 2) is a key upstream regulator of mTORC1, which can be directly phosphorylated and inactivated by Akt and ERK1/2, resulting in mTORC1 activation. Akt also activates mTORC1 by phosphorylating PRAS40 (AKTS1), an mTORC1 inhibitor, hence causing its dissociation from raptor (an mTORC1 subunit). TNFα, Wnt signalling, DNA damage and amino-acids are also known to regulate mTORC1 activity. mTORC1 promotes protein synthesis through phosphorylation of eukaryotic translation initiation factor 4E-binding protein 1 (4E-BP1) and S6 kinase 1 (S6K1), amongst others. mTORC1 also controls the synthesis of lipids necessary for membrane biogenesis and negatively regulates autophagy and lysosomes biogenesis 15 .
One of the most interesting aspects of the Hsu et al. 13 data is that it uses two mTOR inhibitors that are known to have a different functional effect on the kinase. Indeed, the authors used Torin1, which is an ATP-competitive mTOR kinase domain inhibitor, in parallel with rapamycin, which is an allosteric mTORC1 inhibitor. This is particularly relevant since feedback activation of the PI3K-Akt pathway upon mTORC1 inhibition is known to occur, and to be a potential problem for uses of mTOR inhibitor as anti-proliferation agents e.g. in cancer. However relevant, these problems are poorly characterised due to the complex relationships between the activities of the two complexes and their inhibitors. We therefore decided to model the effect of the two inhibitors (rapamycin and Torin1) independently, to see if we could learn something in particular about the effects of mTOR inhibition as part of mTORC1. The resulting networks (6 combined parallel optimisations, parameters as above, displayed with 20% tolerance around the best incoming edges) are displayed on Supplementary Fig. 17.
As we can see on Supplementary Fig. 17, a number of elements are common with our own analyses of mTOR inhibition, such as expected effects on KS6B1, RPTOR, AKT1, and the effect on CDK2 discussed above. More interestingly, a number of elements can only be found under Torin1 inhibition of mTOR, and are therefore presumably consequences of the inhibition of mTORC2 activity. These are: the decreased phosphorylation of AKTS1-T246 via AKT1 (as mentioned above, mTORC2 normally activates Akt which then phosphorylates AKTS1 leading to activation of mTORC1), the decreased phosphorylation of the SGK1 substrate NDRG1-S330 (SGK1 being a canonical downstream effector of mTORC2 activation), and the decreased phosphorylation of AKT3-S472. Of course, a number of elements are also common between the rapamycin and Torin1 treatments, which is expected since both inhibit mTORC1 (and feedbacks between mTORC1 and 2 are complex). Most of these common elements can be directly related to the inhibition of mTORC1. They are: the expected decreased phosphorylation of KS6B1-S447 and RPTOR-S863, the relationship between mTORC1 and DNA damage (as evidenced by the decreased phosphorylation of NPM-S70 and P53-S315), and the decreased phosphorylation of the translation inhibitor PDCD4. More interestingly, there are also a series of elements that are only observed under rapamycin inhibition of mTORC1. These are mostly increases in phosphorylation and were missed by Hsu et al. because increase in phosphorylation are difficult to apprehend when looking at direct effects of kinase inhibition (inhibition is expected to produce decreased phosphorylations), rather than putting phosphorylation perturbations in their network context. These are: the increased phosphorylation of AKT3-S472, of the cell cycle regulator APC1-S688, of the replication associated DNA ligase DNLI1-T233, and of the proteins ELAV1 and F122A. The results in Supplementary Fig. 17B provide evidence for the activation of the MAP kinases cascade and proliferative pathways upon mTORC1, which is not observed under combined mTORC1 and mTORC2 inhibition. This analysis shows how PHONEMeS can be used to analyse the effects of inhibitors with partially overlapping target spectrum. The analysis both recapitulates expected relationships between the inhibitors and their expected targets (clearly placing the known downstream effects of these), and proposes hypotheses for the known phenotypic effects of these inhibitions that have very important consequences for a potential pharmacological use. More generally, this analysis demonstrates the flexibility of our method is flexible and shows that it can in principle be applied to and provide useful insights from any type of untargeted mass spectrometry phosphoproteomics data collected before and after perturbation.