Boolean ErbB network reconstructions and perturbation simulations reveal individual drug response in different breast cancer cell lines
 Silvia Von der Heyde^{1},
 Christian Bender^{2},
 Frauke Henjes^{3},
 Johanna Sonntag^{4},
 Ulrike Korf^{4} and
 Tim Beißbarth^{1}Email author
DOI: 10.1186/17520509875
© von der Heyde et al.; licensee BioMed Central Ltd. 2014
Received: 20 December 2013
Accepted: 10 June 2014
Published: 25 June 2014
Abstract
Background
Despite promising progress in targeted breast cancer therapy, drug resistance remains challenging. The monoclonal antibody drugs trastuzumab and pertuzumab as well as the small molecule inhibitor erlotinib were designed to prevent ErbB2 and ErbB1 receptor induced deregulated protein signalling, contributing to tumour progression. The oncogenic potential of ErbB receptors unfolds in case of overexpression or mutations. Dimerisation with other receptors allows to bypass pathway blockades. Our intention is to reconstruct the ErbB network to reveal resistance mechanisms. We used longitudinal proteomic data of ErbB receptors and downstream targets in the ErbB2 amplified breast cancer cell lines BT474, SKBR3 and HCC1954 treated with erlotinib, trastuzumab or pertuzumab, alone or combined, up to 60 minutes and 30 hours, respectively. In a Boolean modelling approach, signalling networks were reconstructed based on these data in a cell line and time course specific manner, including prior literature knowledge. Finally, we simulated network response to inhibitor combinations to detect signalling nodes reflecting growth inhibition.
Results
The networks pointed to cell line specific activation patterns of the MAPK and PI3K pathway. In BT474, the PI3K signal route was favoured, while in SKBR3, novel edges highlighted MAPK signalling. In HCC1954, the inferred edges stimulated both pathways. For example, we uncovered feedback loops amplifying PI3K signalling, in line with the known trastuzumab resistance of this cell line. In the perturbation simulations on the shortterm networks, we analysed ERK1/2, AKT and p70S6K. The results indicated a pathway specific drug response, driven by the type of growth factor stimulus. HCC1954 revealed an edgetic type of PIK3CAmutation, contributing to trastuzumab inefficacy. Drug impact on the AKT and ERK1/2 signalling axes is mirrored by effects on RB and RPS6, relating to phenotypic events like cell growth or proliferation. Therefore, we additionally analysed RB and RPS6 in the longterm networks.
Conclusions
We derived protein interaction models for three breast cancer cell lines. Changes compared to the common reference network hint towards individual characteristics and potential drug resistance mechanisms. Simulation of perturbations were consistent with the experimental data, confirming our combined reverse and forward engineering approach as valuable for drug discovery and personalised medicine.
Keywords
ErbB RPPA Network reconstruction Boolean model Breast cancer cell line Drug resistanceBackground
Longitudinal time course data are the basis for modelling signalling networks in a holistic systems biology approach in order to uncover mechanisms of signal transduction dynamics [1, 2]. Network models provide novel insight [3, 4] and allow us to perform efficiently simulations to predict systems behaviour or evaluate certain hypotheses [5]. Furthermore, combining perturbation experiments with the measurements of system dynamics seems to be even more efficient than time series data on their own [6–8]. Knockouts or stimuli as directed perturbations support the systematic identification of regulatory relationships.
Quantitative models, based on differential equations, require explicit knowledge on the kinetics of the system of interest [9–12]. In contrast, the qualitative Boolean abstraction considers the components’ states as binary variables, being either active (1) or passive (0), but nevertheless encompasses the essential functionality [13, 14]. Wang et al. stressed, that Boolean models have already been successfully applied in reverse engineering of proteomic signalling networks, and their reduced complexity is considered to be especially advantageous for largescale systems [15]. To avoid the drawbacks of purely data or literaturedriven algorithms regarding completeness, generalisation or interpretability, combined approaches become more and more prominent in the area of network reconstruction [6, 16, 17]. Some reverse engineering approaches, like ddepn[6] or CellNOptR[18], ideally join perturbed time course input data and literature prior knowledge in network reconstruction, while preserving the simplicity of Boolean logic at the same time. Forward engineering methods allow subsequent analysis of the stable states of the reconstructed system. Hence, this may allow to deduce possible longterm behaviour of components activity under perturbations. Such approaches are integrated and freely available in the open source Python software package BooleanNet[19] or in the R[20] package BoolNet[21], for example. As reviewed by Samaga and Klampt [22], several software tools can be applied for the dynamic modeling of logical signal transduction networks. Among others, they exemplarily mentioned GINsim[23], SQUAD[24], BooleanNet[19], ChemChains[25], Odefy[26], and BoolNet[21].
Here we focus on protein signalling networks in breast cancer, representing the most common cancer type among women [27]. Breast cancer, as a heterogeneous disease, can be divided into subgroups, which differ in cellular properties as well as in prognosis. This requires individual therapy approaches, which are in the focus of current research and have partially already been realised.
Here we are interested in the ‘HER2positive’ subtype of breast cancer, overexpressing the human epidermal growth factor receptor 2 (HER2, also termed ErbB2). ErbB2 is a receptor tyrosine kinase (RTK) and member of the epidermal growth factor (EGF) receptor family, consisting of three further RTKs, namely ErbB1, ErbB3 and ErbB4. These receptors cooperatively function as homo or heterodimers after activation via growth factors like EGF for ErbB1 or heregulin (HRG) for ErbB3 [28]. This initialises signalling cascades, pathologically contributing to tumourigenesis and tumour progression. Interestingly, different dimer formations induce different signalling pathways, like PI3K and MAPK, also with differing signalling strengths [29]. The role of the orphan receptor ErbB2 in dysregulation of the ErbB network is of major interest, due to its overexpression in 1020% of breast tumours, diagnosed as HER2positive. Furthermore, its role as favoured dimerisation partner independent on ligandactivation implies oncogenic potential [30–32]. The therapeutic antibodies trastuzumab and pertuzumab have especially been designed to target ErbB2 [33].
However, frequently occurring therapy resistance reduces the efficiency of targeted therapeutics [34–36]. This resistance is often associated with deregulated pathway activity [37, 38] or bypasses via other RTKs, especially ErbB family members [39]. Mainly ErbB1 expression has been anticipated as molecular cause to overcome impact of ErbB2 targeting drugs. Smallmolecule inhibitors such as erlotinib are already in use against nonsmall cell lung cancer [40] and pancreatic cancer [41].
Here we aim as a first step at the identification of individual drug response patterns and insights into drug resistance in HER2positive breast cancer. ErbB2 amplified cell lines were therefore subjected to short and longterm drug treatment with erlotinib, pertuzumab and trastuzumab, alone or in combinations. Samples were analysed by reversephase protein arrays (RPPA) [42]. We were interested in synergistic benefits of combining erlotinib, pertuzumab or trastuzumab in ErbB1 expressing, ErbB2 amplified tumours with differing resistance phenotypes. Therefore three representative breast cancer cell lines were selected as model systems, namely BT474, SKBR3 and HCC1954, of which the latter is known to be trastuzumab resistant due to a PIK3CA mutation, while BT474 exhibits wild type behaviour [43]. The SKBR3 cell line is supposed to be pertuzumab resistant [44].
ErbB dimers predominantly activate the MAPK and PI3K pathway [29]. Therefore, we concentrated on the involved key regulators in fast downstream signalling. Among those were ERK1/2 and AKT, and also p70S6K, which is upstream influenced by both of the signalling axes. Phosphorylation of RPS6 and RB was used as longterm indicator for proliferation, cell cycle or tumour progression [28]. Prior literature knowledge on ErbB signalling was used as input for protein network reconstruction per cell line via ddepn. Beyond that, we inferred combined therapies that target ErbB family members, customised to the topology of the different subtypes. BoolNet was applied to compute stable cycles of protein activity states, socalled attractors, incorporating all possible treatment combinations. This way, optimal drug treatment to deactivate oncogenic proteins was identified.
Methods
Data
Protein abundance and phosphorylation measurements in BT474, SKBR3 and HCC1954 cells were carried out as described by Henjes et al. [28]. In principle, the RPPA protein array technology works as follows. Minimal amounts (1 nl volume) of cell lysate are spotted along with a serial dilution of control samples on nitrocellulosecoated glass slides using a printing robot (Aushon 2470 arrayer). Samples are organised as ordered subarrays so that they are addressable during the data analysis procedure, and a single slide can accommodate one or more subarrays. Each subarray is analysed using a highly specific detection antibody to measure the abundance of a certain protein or its phosphorylation rate. For each spot, the ratio of bound detection antibody is visualised using secondary antibodies labelled with near infrared (NIR) fluorescent dyes. Slides are scanned using the Odyssey scanner (LiCor Biosciences). Spot intensities are determined using a microarray image analysis software (GenePix).
Apart from the quantitative character, another advantage of the technology is the handling of large sample sets which protein abundance can be detected simultaneously in a high throughput fashion. 20200 identical slides can be produced in parallel in a single print run.
In order to normalise the data spotwise for deviant total protein concentrations due to spotting variance, staining with Fast Green FCF dye was employed [42]. Therefore, one slide was stained with the dye to determine the total protein content of each lysate spot and corresponding signal intensity correction factors. The spots on the remaining slides were divided by these correction factors and afterwards multiplied by the median value to scale the data back to the native range.
The RPPA data used here include data presented in Henjes et al. [28]. Additionally, further targets have been measured and were used for network reconstruction. The complete data set has been submitted to the Gene Expression Omnibus (GEO) with accession number GSE50109.
Shortterm measurements
In the shortterm measurements, trastuzumab, pertuzumab and erlotinib were added to the cells in starvation medium one hour before stimulation with the growth factors EGF and HRG. All possible 24 combinations of drugs and stimuli were measured. Application of the stimuli was defined as time point zero in the measurements. The growth factors were chosen to activate explicitly the MAPK and PI3K pathway. Lysate preparation was performed at ten time points, namely after 0, 4, 8, 12, 16, 20, 30, 40, 50 and 60 minutes. The drug treatment experiments comprised three biological replicates, whereas the inhibitorfree experiments incorporated five biological replicates. The experiments for the SKBR3 cell line comprised only two biological replicates of HRG stimulated cells under the triple drug combination. Each biological replicate was spotted in triplicate on the RPPA slides. To obtain shortterm signal intensities, eleven antibodies for specific phosphorylation sites were selected according to quality checks, including inspection of corresponding dilution series and comparison to signals arising from secondary antibodies only. The chosen target proteins and respective antibodies are listed in Additional file 1.
Longterm measurements
For longterm measurements, no explicit ligand stimulation was performed. Instead, cells were incubated in full growth medium for 24 hours prior to adding the three mentioned therapeutics in double combinations or as triplet. Single drug treatment was just conducted with erlotinib. Full growth medium was used to avoid confounding effects of nutrient deficiency. Protein abundance was also quantified without any drug application. The measuring points included 0, 1, 2, 4, 6, 8, 12, 18, 24 and 30 hours with three biological and technical replicates each. At time point 18, only two biological replicates were available. Additional file 1 displays the 21 targets of interest for longterm signalling.
Statistical inference of drug effects
To determine, whether a specific drug treatment revealed an inhibiting effect on the signal intensities of the proteins, we applied the following method. Firstly, for each protein and (combinatorial) drug treatment we linearly modelled the signal intensities as depending on the factors time and group, i.e. no drug treatment versus drug treatment. If the interaction of both factors showed a significant (pvalue < 0.05) influence on the signal intensity, we further applied a Wilcoxon rank sum test for the measurements at time point 60 minutes for the shortterm data, or at time point 30 hours for the longterm data. Thereby, we tested for significantly (pvalue < 0.05) smaller intensity values in the drug treated group. The drug treatments with a significant test result were considered as efficient inhibitors. The therapeutic (combination) with the smallest pvalue was defined as the optimal one.
Literature prior knowledge
We manually determined two reference networks, i.e. one for each time course, as initial joint hypotheses for all of the three breast cancer cell lines. Because emphasis was put on phosphoproteomic signalling, this was mainly based on PhosphoSitePlus^{®;}[45]. Several publications confirm these assumptions, as depicted in Additional file 2.
Network reconstruction
For Boolean network reconstruction, we chose the method of dynamic deterministic effects propagation networks (DDEPN) [6]. This method was particularly tailored to perturbed longitudinal protein phosphorylation data. It is based on the DEPN approach [46], which stands for deterministic effects propagation networks. The determinism is related to the way of perturbation effect propagation in the networks from parent to child nodes, implying transitively closed graphs. The dynamic version of Bender et al. [6, 47] differs with respect to the integration of perturbed time course measurements. While the DEPN approach requires many perturbations, like knockdowns, but only few time points, which are regarded as independent measurements, ddepn is designed for longer time series without the necessity of many or all network nodes being perturbed. The latter situation, i.e. few perturbations by drug interventions, reflected the design of the RPPA experiments under consideration here, hence leading to the application of ddepn. Most network reconstruction algorithms have been designed for gene expression data from microarray measurements [7], which differ from (phospho)protein data regarding the amount of involved network nodes. Many current methods are tailored to the inference of gene regulatory networks based on static measurements at one time point, reflecting the steady state of the system under consideration [48]. The longitudinal time course data used here require a suitable method, as provided by ddepn. The method of Bender et al. was shown to outperform two dynamical Bayesian network approaches, and to be capable of inferring known signalling cascades in the ErbB pathway [47]. A further advantage was the public availability of ddepn as an R[20] package.
The reconstruction procedure is depicted in Additional file 3, and the core elements are described according to [6, 47] in the following. The protein interaction networks are modelled as directed, possibly cyclic, graphs, with nodes V={v _{ i }:i∈1,…,N} representing proteins and edges representing interactions. Also the external perturbations, i.e. the drugs and growth factors in our case, are modelled as nodes. The edge types can be either activating or inhibiting, denoted by 1 and 1, respectively, in the adjacency matrix Φ=V×V→{0,1,−1} of the network. An entry of zero indicates no edge between two nodes. So each edge incorporates a pair of nodes {ϕ _{ i j }:i,j∈1,…,N}. The measurement data, which form the basis for the reconstruction, are stored in a matrix D={d _{ i t r }:i∈1,…,N,t∈1,…,T,r∈1,…,R}, considering T time points and R replicates.
The distribution parameters are for each protein estimated as empirical mean and standard deviation of all measurements for the considered protein in the corresponding class, yielding the parameter matrix $\widehat{\Theta}=\left\{{\widehat{\theta}}_{i0},{\widehat{\theta}}_{i1}\right\}=\left\{({\widehat{\mu}}_{i0},{\widehat{\sigma}}_{i0}),({\widehat{\mu}}_{i1},{\widehat{\sigma}}_{i1})\right\}\phantom{\rule{1em}{0ex}}\forall i\in 1,\dots ,N$.
The prior probability distribution P(Φ) includes penalisation of differences between the network structure Φ and a userdefined prior belief B=V×V→[ −1,1], where the absolute value correlates with the confidence in an edge. Here we chose B=V×V→{0,1,−1}, assuming in advance specific activating, inhibiting or missing edges with maximum confidence. We made use of the Laplace prior model (laplaceinhib), accounting for both edge types, i.e. activation and inhibition. The prior belief for an edge is defined as $P({\varphi}_{\mathit{\text{ij}}}{b}_{\mathit{\text{ij}}},\lambda ,\gamma )=\frac{1}{2\lambda}{e}^{\frac{{\Delta}_{\mathit{\text{ij}}}}{\lambda}}$, including a weighted difference term Δ _{ i j }=ϕ _{ i j }−b _{ i j }^{ γ } with a weight exponent $\gamma \in {\mathbb{R}}^{+}$. As the edge probabilities are assumed to be independent, the prior belief for a network structure Φ is derived as the product of those, i.e. $P(\Phi B,\lambda ,\gamma )=\prod _{i,j}P({\varphi}_{\mathit{\text{ij}}}{b}_{\mathit{\text{ij}}},\lambda ,\gamma ),\phantom{\rule{1em}{0ex}}i,j\in \left\{1,\dots ,N\right\}$. The individual edge probabilities lie between 0 and $\frac{1}{2\lambda}\phantom{\rule{1em}{0ex}}\forall \lambda ,\gamma \in {\mathbb{R}}^{+}$. The protein interactions corresponding to our chosen prior are displayed in Additional file 2. The prior’s impact strength was emphasised in such a way, that only strongly deviating data influence the network structure, because the ErbB wiring as well as the MAPK and PI3K pathways are well examined in literature. This prioritisation is reflected in the hyperparameter λ set to 0.0001. For the parameter γ we chose one, neglecting extra penalisation of deviation from the prior. These settings should preserve robustness, but at the same time allow enough impact strength of strongly differing data values.
The network inference via inhibMCMC spanned 50,000 iterations with the first 25,000 iterative steps as burnin phase. To ensure convergence, ten parallel MCMC chains were run, each initialised with a starting network. Convergence was validated via Gelman diagnostic [49]. Nine of the initial ten networks were randomly generated, i.e. for the defined nodes activating, inhibiting or no edges were sampled. The remaining network assumed no connections between the nodes. These initial networks were pruned to the following constraints. Firstly, the nodes related to the growth factors and drugs must not have any ingoing edges. Above that, the indegree of all nodes was limited to four. Finally, no selfloops were allowed. To find significantly occurring edges among the independent runs, merging into a consensus network, a Wilcoxon rank sum testing procedure was used. In detail, in each run the amount of sampled activations and inhibitions per edge was counted and divided by the total number of sampled edges. Subsequently the nullhypothesis was tested, whether the means of these ten edgespecific confidence values equal the same for activation and inhibition. In case of not rejecting the nullhypothesis, coming along with an adjusted pvalue exceeding the significance level α=0.05, no edge was assumed. Otherwise, the respective alternative determined the type of interaction. Adjustment for multiple testing followed the method of Benjamini and Hochberg, controlling the false discovery rate [50]. The whole procedure was embedded into a leaveoneout crossvalidation approach. So each of the ten MCMC chains was left out once, and the testing algorithm was applied to the remaining runs. An edge was included in the final consensus network if it occurred in all of the cross validation runs. Finally, to prevent excessive spurious or obsolete connections ascribable to transitivity, as argued by Bo Na Ki et al. [51], newly reconstructed edges were successively added to the prior network according to ddepn significance and fit of resulting attractor states to the observations of Henjes et al. [28].
Perturbation simulations
To figure out which input of drug combination leads to a certain attractor state of the reconstructed network system, the R package BoolNet[21] was applied. The motivation was based on the assumption that attractors, representing cycles of states, comprise the stable states of cell function. In those states networks mostly reside. Hence, they mirror system phenotypes, dependent on the perturbation context. To the best of our knowledge, apart from BoolNet, there are hardly any R packages offering attractor calculations for Boolean networks. This package supports import of networks in form of files containing Boolean formulas. So it could be easily integrated in our workflow as subsequent analysis step after network reconstruction.
We used its functionality to identify attractors in a synchronous and an asynchronous way. The resulting attractors were steadystate attractors. These consist of only one state, in which all transitions from this state result. These attractors are identical for synchronous and asynchronous updates. We focused on the steadystates, as these should reflect the homoeostatic system state of the cell lines. Intermediate transition states would be interesting as well, but due to the large amount of the involved targets, it would have been too complex to analyse those here in detail.
The search started from predefined initial states of the network nodes. The drug and growth factor nodes were fixed to specific values, reflecting the conducted experiment to be simulated. For shortterm signalling, perturbations included all possible combinations of the therapeutics under the combined stimulus of EGF and HRG. Although the data of separate stimulation with EGF and HRG was used for network reconstructions, here we focused on the combined treatment, representing a more natural tumour environment than a single growth factor alone. Two possible binary states, i.e. active (1) or passive (0), to the power of three different drugs led to eight possible combinations. These were used as fixed input conditions, as the effect was assumed to be continuously valid. Analogously, the growth factors were permanently fixed to one. The remaining protein activity start states were initialised with zero. These components were flexible towards updates. In the longterm measurements, no growth factors were involved but full growth medium. This was defined as one stimulating input S, initially activating the ErbB receptors. This also led to eight fixed input combinations.
BoolNet expects network representation in form of logical interaction rules as input. In contrast, ddepn delivers network reconstruction output in terms of adjacency matrices. Therefore, we incorporated an interface function into the ddepn package, called adjacencyMatrix_to_logicalRules. In detail, the loadNetwork function of BoolNet requires a file containing rowwise logical activation rules of each network node. Each row looks like ‘target node, (activator_1  activator_2) & !(inhibitor_1  inhibitor_2)’, here exemplary for a node with two ingoing activating and inhibiting edges each. The logical OR operator is encoded by ‘ ’, the logical AND is encoded by ‘&’, and logical negation is represented by ‘!’. Accordingly, all of the A inferred activating nodes V _{+}={v _{ a }:a∈1,…,A, A<N} of a target node v _{ j }, represented by an adjacency matrix entry ϕ _{ a j }=1, and v _{ j } itself were connected via OR operators. This ensured that at least one of the activators or the target protein itself had to be active to activate the target node. Analogously, the I inhibiting nodes V _{−}={v _{ i }:i∈1,…,I, I<N} with ϕ _{ i j }=−1 were connected via OR operators. A logical negation operator was attached to ensure that the activity of one of the nodes v _{ i } would result in an inactive node v _{ j }. Both sets of activators and negated inhibitors were then connected via a logical AND operator. After conversion of the adjacency matrices to logical rules, those were implemented in BoolNet into a computational model, to perform perturbation simulations per cell line and time course as well as subsequent analyses of the resulting attractor states.
Results and discussion
Shortterm signalling network reconstruction
HCC1954 is driven by the PI3K as well as the MAPK pathway
In HCC1954, the new edges contributed to both, PI3K and MAPK, signalling. The interaction ErbB1 →ErbB2 reflected a dominant role of heterodimerisation of both receptors, as described by Henjes et al. [28]. The fact that it was specifically inferred for HCC1954, pointed to hyperactive ErbB1/2 heterodimers here. These are known to trigger the MAPK but also, to a lesser extent, the PI3K pathway. The link PDK1 →MEK1/2, supported by Sato et al. [55], stressed crosstalk between these pathways, placing PDK1 into a key position in the PI3K pathway, and MEK1/2 in the MAPK pathway, respectively. Two of the new edges in HCC1954, PDK1 →ErbB2 and p70S6K →AKT, contributed to feedback loops, which were not present in the other two cell lines. Such a topological network element could stabilise the known trastuzumab resistance by boosting the oncogenic effect of ErbB2 and the mutant hyperactive PI3K pathway. Evidence for the feedback mechanism involving PDK1 was provided by Maurer et al. [56] and Tseng et al. [57]. Vega et al. noted an indirect activation of AKT by p70S6K via mTOR [58].
BT474 is driven by the PI3K pathway, while SKBR3 is driven by the MAPK pathway
Comparably to HCC1954, in BT474 an edge indicating hyperactive heterodimers was found, namely ErbB3 →ErbB2, here interestingly with a strong impact on AKT [28]. BT474 is known to contain a rare type of PIK3CA mutation [43]. Pathway crosstalk was also observed in BT474, but here MEK1/2 activated PDK1, and not vice versa like in HCC1954. This edge was supported by Frödin et al. [59], underlining dominant PI3K signalling in this cell line.
The newly detected interactions in SKBR3 started from ErbB3 and PDK1, and both activated ERK1/2. This reflected a dominant MAPK pathway, in which ErbB3 →ERK1/2 was interpretable as indirect stimulation of ERK1/2 via MEK1/2, activated by ErbB2/3 dimers [55].
Perturbation simulations on shortterm networks
Attractor states of shortterm perturbation simulations
BT474  HCC1954  SKBR3  

Simulation  AKT  ERK1/2  p70S6K  AKT  ERK1/2  p70S6K  AKT  ERK1/2  p70S6K  
A  E  A  E  A  E  A  E  A  E  A  E  A  E  A  E  A  E  
X  1  1  1  1  1  1  1  1  1  
E  0  1  1  0  1  0  1  0  1  0  1  0  0  1  
P  1  1  1  0  1  0  1  0  1  0  1  0  1  0  1  
T  1  1  1  1  1  1  1  1  1  
E, P  1  1  0  1  0  1  0  0  1  0  0  0  0  1  
E, T  0  1  1  0  1  1  0  1  1  0  1  0  0  1  
P, T  1  1  1  0  1  0  0  1  0  0  0  1  
E, P, T  1  1  0  1  0  1  0  0  1  0  0  0 
Optimal drug treatment in shortterm signalling
Cell line  AKT  ERK1/2  p70S6K 

BT474    PE  PE 
HCC1954  PT E  E  P E 
SKBR3  PTE  E  PTE 
Inhibition of PI3K signalling requires drug combinations
In SKBR3, the triple drug combination was most effective in inhibiting AKT (Figure 4, Table 2). In BT474, pertuzumab combined with erlotinib was most efficient, but AKT signalling was not fully suppressed as in SKBR3 (Figure 5). Statistically, we did not infer any significant positive drug effect in this cell line. Obviously, erlotinib in synergistic combination with at least pertuzumab was needed to block the ErbB2 receptor and its heterodimerisation, mainly with ErbB1, but also ErbB3. The HRG activated ErbB2/3 heterodimers and PI3K pathway in BT474, as revealed by the network reconstructions, might have prevented a potent drug efficacy.
Interestingly, BT474 and SKBR3 required pertuzumab. This drug was especially designed to prevent heterodimerisation with ErbB2. The stimuli EGF and HRG together activate PI3K signalling by ErbB2/3, ErbB1/2 and ErbB1/3 dimers (Figure 2). The need for pertuzumab combined with erlotinib indicated an important role of ErbB1/2 dimers. This was supported by the fact, that in HCC1954 with dominant heterodimers of this type, as revealed by network reconstructions, none of the drugs was likewise efficient in inhibiting AKT (Figure 6). However, the optimal effect was revealed for the triple drug combination (Table 2). The simulations suggested pertuzumab alone or a combination of both monoclonal antibodies (Table 1). It has to be kept in mind, that the attractor states resembled a longterm steady state, which can differ from observations up to 60 minutes.
The perturbation simulations in BT474 did not lead to inactive AKT upon combined pertuzumab and erlotinib treatment. Instead, erlotinib alone or combined with trastuzumab was efficient (Table 1). Nevertheless, this supported the need for the small molecule inhibitor and a monoclonal antibody to suppress ErbB2 induced PI3K signalling. In SKBR3, the attractor states confirmed the described optimal drug treatment to deactivate AKT. Trastuzumab, when applied alone, was the only treatment without a positive effect in the simulations (Table 1).
Inhibition of MAPK signalling requires erlotinib
Signalling through the MAPK pathway, represented by ERK1/2 activation, was efficiently inhibited by erlotinib alone in both, HCC1954 (Figure 6, Table 2) and SKBR3 (Figure 4, Table 2), cell lines. EGF activates the MAPK pathway via ErbB1 homodimers and ErbB1/2 heterodimers (Figure 2). Both are prevented by ErbB1 inhibition via erlotinib, which was especially designed to target this receptor.
In BT474, pertuzumab plus erlotinib was required (Figure 5, Table 2). This was analogous to the situation in PI3K signalling.
HRG activates the MAPK pathway via ErbB2/3 heterodimers (Figure 2). Obviously, BT474 needed the addition of the monoclonal antibody due to dominant ErbB2/3 formation and activity. On the contrary, the other two cell lines just needed erlotinib alone. Here, in addition to the ErbB1 dimers, the ligandindependent ErbB2 homodimers might have driven ERK1/2 activation and could be inhibited by the small molecule inhibitor. Efficacy of erlotinib towards ErbB2 dimers was previously mentioned by Schaefer et al. [61].
In BT474, the simulations resulted in active ERK1/2 states, resisting drug treatment (Table 1). In HCC1954 and SKBR3, the positive effect of erlotinib was supported by the simulations. The attractor states were additionally inactive for all other (combinatorial) drug treatments, but not trastuzumab alone.
p70S6K is influenced by both, PI3K and MAPK, pathways
The target p70S6K is upstream influenced by the PI3K as well as the MAPK pathway (Figure 2). Hence, p70S6K merges both pathways, leading to activation of RPS6 [60].
The three cell lines showed different pathway preferences. BT474 required the combination of pertuzumab and erlotinib to suppress p70S6K (Figure 5). On the contrary, in SKBR3 the triple drug combination was shown to be optimal (Table 2). Obviously, the effect was driven by erlotinib (Figure 4), which was supported by the attractor states of p70S6K (Table 1). This resembled the drug response of ERK1/2 and reflected a stronger influence by the MAPK pathway. In HCC1954, deactivation of p70S6K was reached via application of erlotinib combined with pertuzumab (Table 2). The treatment with erlotinib alone had a similar effect (Figure 6), while the simulations just confirmed a positive effect of pertuzumab (Table 1). Thus, this cell line seemed to be influenced by both, PI3K and MAPK, pathways.
These results were in line with the newly inferred edges in the network reconstructions. They pointed to a strong influence of PI3K in BT474 in contrast to a dominant MAPK pathway in SKBR3. HCC1954 was influenced by both pathways to a similar extent.
To follow up on the hypothesis that different pathways contribute to a different extent in individual cell lines, we tested correlation between the p70S6K time course and the ones of AKT and ERK1/2, respectively. In BT474, p70S6K correlated positively with AKT (pvalue 0.01, Kendall’s τ estimate 0.64). In HCC1954, p70S6K correlated positively with both, AKT (pvalue <2.22·10^{−16}, Kendall’s τ estimate 0.69) and ERK1/2 (pvalue <2.22·10^{−16}, Kendall’s τ estimate 0.87). In SKBR3, p70S6K also correlated positively with both, AKT (pvalue 0.05, Kendall’s τ estimate 0.51) and ERK1/2 (pvalue 0.02, Kendall’s τ estimate 0.6), with a stronger tendency towards MAPK signalling. The correlation was not as convincing as in the other two cell lines. One could speculate, that the dominance of the MAPK pathway in SKBR3 cells was not as strong as the dominance of the PI3K pathway in BT474. This was supported by the reconstructed networks. They revealed downstream effects of MAPK signalling in SKBR3, while they revealed hyperactive ErbB2/3 dimers in BT474. The dimers drive PI3K already at the receptor layer, and especially ErbB2/3 dimers are regarded as the most potent heterodimer [29].
Drug resistance in HCC1954 regarding the PI3K pathway
In HCC1954, the inferred optimal treatment against AKT signalling with the triple drug combination was not convincing (Figure 6). Analogously, Henjes et al. did not monitor any positive drug effect on AKT under EGF application alone [28]. However, the simulations suggested pertuzumab alone or a combination of both monoclonal antibodies to inhibit AKT phosphorylation. In principle, divergence of simulations from experimental observations can be expected, as the simulated steady state of the system does not necessarily have to be reached after the measured period of time. Anyhow, the apparent resistance here pointed to a hyperactive PI3K pathway which was explainable by the newly inferred HCC1954 edges described in the previous subsection. They represented feedback loops, hyperactive ErbB1/2 heterodimers and pathway crosstalk. On the contrary, in SKBR3, the triple drug combination worked well, as described before. The simulations even predicted efficacy of every other drug (combination) apart from trastuzumab alone. The drug efficacy towards AKT in this cell line could be explained by the fact that the two reconstructed interactions in SKBR3 mainly promoted the MAPK instead of the PI3K pathway.
The regulation of AKT activity under drug influence, highly diverging in HCC1954 and SKBR3, attracted our attention. Therefore we intended testing for edgetic mutations, as discussed by Zhong et al. [62], leading to AKT gainoffunction in HCC1954. Such mutations, perturbing not a node but an edge of a network, are speculated to have deeper impact on phenotypic manifestation of a disease. In detail, we removed each of the AKT stimulating edges outgoing from p70S6K, PDK1, mTOR and ErbB3, alone or in all possible eleven combinations. We then computed the attractor states for the modified networks in HCC1954.
Removal of the connections of mTOR, PDK1 and ErbB3 alone or combined had no influence on improving drug effects, i.e. AKT just got inactive under pertuzumab treatment. Involvement of p70S6K →AKT in the withdrawal process led to much better results. Removed alone or in double combinations with the aforementioned edges, as well as in the two triple combinations containing mTOR, AKT was deactivated under all drug treatments, but not yet trastuzumab alone. Finally, simultaneous removal of the outgoing connections from p70S6K, ErbB3 and PDK1 with or without mTOR, turned out as the only combination enabling potency of all possible drug combinations, including trastuzumab alone. This hinted at a less strong impact of mTOR on AKT here, but indicated synergistic drug resistance potential of p70S6K, ErbB3 and PDK1, also due to the newly inferred edges.
Longterm signalling network reconstruction
The reconstructed longterm signalling networks per cell line are displayed in the Additional file 6. Additional file 5 lists the equivalent Boolean logical interaction rules. Compared to the prior network, most of the newly inferred edges were individual for each cell line, but HCC1954 shared ErbB1 →ERK1/2 with SKBR3, for example. This seemed to be an indirect edge via cRAF, as represented in the prior network. Besides activating connections, also inhibiting ones and edge deletions occurred. For HCC1954, ten new interactions were reconstructed, while two were deleted. In BT474, nine new links were added, and one edge was deleted. In SKBR3, we inferred 20 new connections and one deletion, namely the removal of p53 activation via p38, bearing oncogenic risk [63, 64].
In contrast to the shortterm networks, new feedback loops were reconstructed in every cell line, not exclusively in HCC1954. In HCC1954, the mutual activation between p53 and RB established such a feedback mechanism. For SKBR3 we even inferred two edges, each forming feedback loops. Contrary to HCC1954, p53 inhibited RB. The second loop connection was inhibition of ErbB3 by AKT, pointing to a negative feedback against PI3K signalling [65–67].
In HCC1954, the newly inferred edges Cyclin B1 →AKT and ErbB3 →ErbB1 contributed to PI3K signalling, of which the latter was explainable as heterodimers. The newly inferred edge cJUN →ErbB1 in HCC1954 also indicated raised activity of ErbB1. Interestingly, in SKBR3 we conducted an inhibiting edge from Cyclin B1 to AKT but instead an activating one to ERK1/2, contributing to MAPK signalling, which was also stated by Abrieu et al. [68]. Another new edge in HCC1954 involved a cell cycle player, i.e. activation of Cyclin D1 by p70S6K [69]. Accordingly, we inferred RPS6 →Cyclin D1 in BT474, with RPS6 as downstream target of p70S6K. In SKBR3, the edge p70S6K →Cyclin B1 was reconstructed. A further interesting new activating edge in HCC1954 led from RB to TSC2, while we inferred a reversed inhibition in SKBR3. Searle et al. discussed targeting RB deficient cancers by deactivating TSC2 [70].
Two novel interactions in BT474 activated Cyclin B1, arising from ErbB1 and ErbB3, respectively, which meant that mitosis was driven by ErbB1/3 dimers in this cell line. This indicated a hyperactive PI3K pathway, as revealed in the shortterm case.
In SKBR3, we reconstructed an outgoing edge from the artificial network stimulus S, representing full growth medium, activating AKT. This could be explained as strong activation of AKT, driving PI3K signalling in this cell line. The new edges ErbB2 →TSC2 and ErbB3 →PRAS had to be interpreted as indirect effects, too. They pointed to activity of ErbB2/3 dimers, feeding into both, MAPK and PI3K, pathways. The edge ErbB2 →TSC2 could imply an oncogenic role of TSC2. Liu et al. discussed a context dependent functionality of TSC2 [71].
Perturbation simulations on longterm networks
Attractor states of longterm perturbation simulations
BT474  HCC1954  SKBR3  

Simulation  RPS6  RB  RPS6  RB  RPS6  RB  
A  E  A  E  A  E  A  E  A  E  A  E  
X  1  0  1  1  0  1  1  0  1  
E  0  1  0  0  1  0  0  1  0  
P  0    0    0    0    0    0   
T  1    0    1    0    1    0   
E, P  0  1  0  0  1  0  0  1  0  
E, T  0  1  0  0  1  0  0  0  
P, T  0  0  0  1  0  1  0  0  
E, P, T  0  1  0  0  1  0  1  0  0 
The control measurements without any drug treatment should result in activation of ErbB members and dimerisation events, promoting cell growth and proliferation. In fact, this was expressed as reasonable activation of AKT, ERK1/2 and RPS6 in all cell lines, which held for simulations as well as experimental observations. In contrast, the attractor states of RB were inactive in all cell lines (Table 3). Actually, a continuously rising stimulation effect over 30 hours was not observed for HCC1954 and SKBR3 by Henjes et al. [28] either.
Optimal drug treatment in longterm signalling
Cell line  AKT  ERK1/2  RB  RPS6 

BT474    TE  E  TP 
HCC1954  E  TE  E   
SKBR3  PTE  TE  TE  TE 
The optimal longterm drug response for AKT and ERK1/2 confirms shortterm observations
This was also the explanation, why we detected erlotinib, but not the combination with trastuzumab, as the optimal treatment in HCC1954 (Table 4). In SKBR3, we inferred the triple drug combination as the optimal one, but the combination of both monoclonal antibodies alone also had a significant effect over time (Figure 7). Hence, like in the shortterm results, a drug combination was required to suppress PI3K signalling, here with an obvious need for trastuzumab. For BT474 and HCC1954, this was supported by the simulation results, in which trastuzumab alone had no effect, but was efficient within drug combinations. In HCC1954, even the best drug response was not convincing (Figure 7), pointing to a dominant PI3K pathway, as revealed in the shortterm analysis. Interestingly, SKBR3 showed a strong activation peak of AKT phosphorylation between 8 and 18 hours (Figure 7), which was just suppressed under combined application of trastuzumab and pertuzumab. We revealed a positive correlation with ERK1/2 (pvalue 0.02, Kendall’s τ estimate 0.6) and RPS6 (pvalue 0.01, Kendall’s τ estimate 0.64). The reconstructed edges S →AKT and ErbB1 →ERK1/2 in SKBR3 indicated strong activation of AKT and ERK1/2. In addition to the prior network, in which AKT and ERK1/2 fed into RPS6 phosphorylation via p70S6K, some of the novel edges pointed to a positive feedback from p70S6K or RPS6 to ERK1/2. The feedback from p70S6K via Cyclin B1, for example, was expressed by the edges p70S6K →Cyclin B1 and Cyclin B1 →ERK1/2. Compared to the shortterm results, indicating a dominant MAPK pathway, this longterm observation indicated strong signalling via both, PI3K and MAPK, pathways in SKBR3.
Quick drug response for RPS6 and delayed response for RB
In HCC1954, it was the combination of both monoclonal antibodies, that failed in deactivating RPS6 (Figure 9), while the simulations predicted trastuzumab alone to fail (Table 3). The graphical observations were similar for RB (Figure 10). The newly inferred edges ErbB3 →ErbB1 and cJUN →ErbB1 in HCC1954 explained the necessity for erlotinib against ErbB1 dimers. The positive impact of erlotinib, the optimal treatment against RB signalling (Table 4), was supported by simulations. However, it did not unfold before 1218 hours, in case of RB as well as RPS6. Regarding RPS6, no significant effect was detected for HCC1954 (Table 4).
According to Henjes et al. [28], in SKBR3 erlotinib and all therapeutic combinations helped to suppress RPS6, which was supported by the simulations (Table 3). As shown in Figure 9, the combination of trastuzumab and erlotinib was the only one, that revealed its continuous inhibiting effect already after one hour. This combined treatment was also statistically inferred as the optimal one (Table 4). The same combination was optimal with respect to RB activity, which was also in line with the simulations. Here, analogously to BT474 and HCC1954, the drug effect did not appear before 18 hours (Figure 10).
As the combination of trastuzumab and erlotinib was efficient in all of the three cell lines against RPS6 as well as RB phosphorylation, we further analysed target correlations under this drug combination to explain the different rapidness of drug responses.
In BT474, RB positively correlated with Cyclin B1 (pvalue 0.02, Kendall’s τ estimate 0.6), while RPS6 positively correlated with ERK1/2 (pvalue 0.02, Kendall’s τ estimate 0.6). Obviously, RPS6 was mainly stimulated by the MAPK pathway, which was efficiently inhibited by the combination of trastuzumab and erlotinib in a fast manner. On the contrary, RB seemed to be influenced by Cyclin B1. The newly reconstructed edges ErbB1 →Cyclin B1 and ErbB3 →Cyclin B1 supported hyperactivity of Cyclin B1, driven by ErbB1/3 heterodimers. In SKBR3, RB negatively correlated with PRAS (pvalue 0.05, Kendall’s τ estimate 0.51) and TSC2 (pvalue 0.03, Kendall’s τ estimate 0.56), while RPS6 positively correlated with AKT (pvalue 0.03, Kendall’s τ estimate 0.56) and ERK1/2 (pvalue 0.02, Kendall’s τ estimate 0.6). Obviously, like in BT474, RPS6 was mainly activated through the MAPK pathway. Interestingly, RB seemed to require inhibition via PRAS or TSC2. The latter was confirmed via one of the novel edges in SKBR3, namely inhibition of RB by TSC2. In addition, PRAS as well as TSC2 seemed to be especially active in this cell line with regard to the new edges ErbB3 →PRAS and ErbB2 →TSC2.
In HCC1954, the drug response was not only delayed for RB, but also for RPS6, which was in line with the positive correlation with RB (pvalue <2.22·10^{−16}, Kendall’s τ estimate 0.73). Like in BT474, Cyclin B1 seemed to be a driving force, since both, RPS6 (pvalue 0.03, Kendall’s τ estimate 0.56) and RB (pvalue <2.22·10^{−16}, Kendall’s τ estimate 0.82) positively correlated with this target. The new edge Cyclin B1 →AKT supported special activation of RPS6 via PI3K signalling, leading to a delayed drug response. Interestingly, we revealed negative correlations, as observed for SKBR3. In HCC1954, RPS6 and RB correlated with BAX (pvalue 0.03, Kendall’s τ estimate 0.56) and FoxO1/3a (pvalue 0.05, Kendall’s τ estimate 0.51), pointing to a delayed inhibition of RPS6 and RB via BAX or FoxO1/3a.
Conclusions
Using a combination of reverse and forward engineering techniques, we focused on deregulated protein interactions in the ErbB network in a Boolean modelling framework. The reconstructed hypothetical networks revealed individual protein interactions contributing to signalling pathway preferences as well as drug resistance via feedback loops, pathway crosstalk or hyperactive heterodimers. While this reverse engineering focused on the network edges, we concentrated in the subsequent forward engineering step on the network nodes. The perturbation simulations for AKT, ERK1/2, RB and RPS6 mainly confirmed our graphical and statistical analyses as well as the observations of Henjes et al. [28] regarding (combinatorial) drug efficacy. However they have to be interpreted as an independent, more prospective investigation, because stable system states do not necessarily have to be reached in temporally limited observations.
In the first step, the combined Boolean modelling approach revealed the mechanisms underlying individual drug response. In the second step, it predicted the network propagation effects on protein activity, and hence the drug response itself.
One major finding is, that different breast cancer phenotypes seem to be driven by specific pathway preferences in the ErbB network. This leads to individual drug response, requiring different therapeutic treatments. The perturbation simulations revealed a more diverse drug response in shortterm than in longterm signalling, which stresses the importance of early intervention at the top level layer of the signalling network.
Another interesting aspect is to combine edge and node perturbations in Boolean network models to reveal edgetic mutations, as we did in the HCC1954 cell line for AKT.
Basic molecular research, embedded in a Boolean modelling framework here, composes a first step to gain insight into individual mechanisms of drug response or resistance mechanisms in breast cancer. Especially, the proteomic signalling interplay directly effects tumour development and represents a promising target in cancer therapy, which has to be understood in more detail in the future.
Abbreviations
 DDEPN:

Dynamic deterministic effects propagation networks
 EGF:

Epidermal growth factor
 EGFR:

Epidermal growth factor receptor
 GEO:

Gene expression omnibus
 HER2:

Human epidermal growth factor receptor 2
 HMM:

Hidden Markov model
 HRG:

Heregulin
 MAPK:

Mitogenactivated protein kinase
 MCMC:

Markov chain Monte Carlo
 NIR:

Near infrared
 PI3K:

Phosphoinositide 3kinase
 RPPA:

Reverse phase protein array
 RTK:

Receptor tyrosine kinase.
Declarations
Acknowledgements
Pertuzumab, trastuzumab and erlotinib were provided by Roche Diagnostics GmbH, Penzberg, Germany. This work was supported by a grant from the German Federal Ministry of Education and Research (BMBF) within the Medical Systems Biology programme BreastSys.
We also acknowledge support by the Open Access Publication Funds of the Göttingen University.
Authors’ Affiliations
References
 Hill SM, Lu Y, Molina J, Heiser LM, Spellman PT, Speed TP, Gray JW, Mills GB, Mukherjee S:Bayesian inference of signaling network topology in a cancer cell line. Bioinformatics (Oxford, England). 2012, 28 (21): 28042810. PMID: 22923301,View ArticleGoogle Scholar
 Park Y, Bader JS:How networks change with time. Bioinformatics (Oxford, England). 2012, 28 (12): 4048. PMID: 22689777,View ArticleGoogle Scholar
 Roukos DH:Trastuzumab and beyond: sequencing cancer genomes and predicting molecular networks. Pharmacogenom J. 2011, 11 (2): 8192. PMID: 20975737,View ArticleGoogle Scholar
 Oda K, Matsuoka Y, Funahashi A, Kitano H:A comprehensive pathway map of epidermal growth factor receptor signaling. Mol Syst Biol. 2005, 1: 2005.0010PMID: 16729045,PubMed CentralView ArticlePubMedGoogle Scholar
 Feiglin A, Hacohen A, Sarusi A, Fisher J, Unger R, Ofran Y:Static network structure can be used to model the phenotypic effects of perturbations in regulatory networks. Bioinformatics (Oxford, England). 2012, 28 (21): 28112818. PMID: 22923292,View ArticleGoogle Scholar
 Bender C, Heyde SV, Henjes F, Wiemann S, Korf U, Beissbarth T:Inferring signalling networks from longitudinal data using sampling based approaches in the rpackage ‘ddepn’. BMC Bioinformatics. 2011, 12: 291PMID: 21771315,PubMed CentralView ArticlePubMedGoogle Scholar
 Penfold CA, BuchananWollaston V, Denby KJ, Wild DL:Nonparametric bayesian inference for perturbed and orthologous gene regulatory networks. Bioinformatics. 2012, 28 (12): 233241. PMID: 22689766,View ArticleGoogle Scholar
 Wagner JP, WolfYadlin A, Sevecka M, Grenier JK, Root DE, Lauffenburger DA, MacBeath G:Receptor tyrosine kinases fall into distinct classes based on their inferred signaling networks. Sci Signaling. 2013, 6 (284): 58View ArticleGoogle Scholar
 Chen WW, Schoeberl B, Jasper PJ, Niepel M, Nielsen UB, Lauffenburger DA, Sorger PK:Inputoutput behavior of ErbB signaling pathways as revealed by a mass action model trained against dynamic data. Mol Syst Biol. 2009, 5: PMID: 19156131Google Scholar
 Hatakeyama M, Kimura S, Naka T, Kawasaki T, Yumoto N, Ichikawa M, Kim JH, Saito K, Saeki M, Shirouzu M, Yokoyama S, Konagaya A:A computational model on the modulation of mitogenactivated protein kinase (MAPK) and akt pathways in heregulininduced ErbB signalling. Biochem J. 2003, 373 (Pt 2): 451463. PMID: 12691603,PubMed CentralView ArticlePubMedGoogle Scholar
 Jones RB, Gordus A, Krall JA, MacBeath G:A quantitative protein interaction network for the ErbB receptors using protein microarrays. Nature. 2006, 439 (7073): 168174. PMID: 16273093,View ArticlePubMedGoogle Scholar
 Schoeberl B, EichlerJonsson C, Gilles ED, Müller G:Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors. Nat Biotechnol. 2002, 20 (4): 370375. PMID: 11923843,View ArticlePubMedGoogle Scholar
 Sahin O, FrÃűhlich H, LÃűbke C, Korf U, Burmester S, Majety M, Mattern J, Schupp I, Chaouiya C, Thieffry D, Poustka A, Wiemann S, Beissbarth T, Arlt D:Modeling ERBB receptorregulated G1/S transition to find novel targets for de novo trastuzumab resistance. BMC Syst Biol. 2009, 3: 1PMID: 19118495,PubMed CentralView ArticlePubMedGoogle Scholar
 Samaga R, SaezRodriguez J, Alexopoulos LG, Sorger PK, Klamt S:The logic of EGFR/ErbB signaling: theoretical properties and analysis of highthroughput data. PLoS Comput Biol. 2009, 5 (8): 1000438PMID: 19662154,View ArticleGoogle Scholar
 Wang RS, Saadatpour A, Albert R:Boolean modeling in systems biology: an overview of methodology and applications. Phys Biol. 2012, 9 (5): 055001View ArticlePubMedGoogle Scholar
 McDermott JE, Wang J, Mitchell H, WebbRobertson BJ, Hafen R, Ramey J, Rodland KD:Challenges in biomarker discovery: Combining expert insights with statistical analysis of complex omics data. Expert Opin Med Diag. 2013, 7 (1): 3751. PMID: 23335946,View ArticleGoogle Scholar
 Eduati F, De Las Rivas J, Di Camillo B, Toffolo G, SaezRodriguez J:Integrating literatureconstrained and datadriven inference of signalling networks. Bioinformatics (Oxford, England). 2012, 28 (18): 23112317. PMID: 22734019,View ArticleGoogle Scholar
 Terfve C, Cokelaer T, Henriques D, Goncalves E, Morris MK, van Iersel M, Lauffenburger DA, SaezRodriguez J, MacNamara A:CellNOptR: a flexible toolkit to train protein signaling networks to data using multiple logic formalisms. BMC Syst Biol. 2012, 6: 133PMID: 23079107,PubMed CentralView ArticlePubMedGoogle Scholar
 Albert I, Thakar J, Li S, Zhang R, Albert R:Boolean network simulations for life scientists. Source Code Biol Med. 2008, 3: 16PMID: 19014577,PubMed CentralView ArticlePubMedGoogle Scholar
 R Core Team: R: A Language and Environment for Statistical Computing. 2012, Vienna, Austria: R Foundation for Statistical Computing, R Foundation for Statistical Computing. ISBN 3900051070. http://www.rproject.org.Google Scholar
 Müssel C, Hopfensitz M, Kestler HA:BoolNet–an r package for generation, reconstruction and analysis of boolean networks. Bioinformatics (Oxford, England). 2010, 26 (10): PMID: 20378558View ArticleGoogle Scholar
 Samaga R, Klamt S:Modeling approaches for qualitative and semiquantitative analysis of cellular signaling networks. Cell Commun Signaling: CCS. 2013, 11 (1): 43PMID: 23803171 PMCID: PMC3698152,PubMed CentralView ArticlePubMedGoogle Scholar
 Gonzalez AG, Naldi A, Sánchez L, Thieffry D, Chaouiya C:GINsim: a software suite for the qualitative modelling, simulation and analysis of regulatory networks. Bio Syst. 2006, 84 (2): 91100. PMID: 16434137,Google Scholar
 Di Cara A, Garg A, De Micheli G, Xenarios I, Mendoza L:Dynamic simulation of regulatory networks using SQUAD. BMC Bioinformatics. 2007, 8: 462PMID: 18039375 PMCID: PMC2238325,PubMed CentralView ArticlePubMedGoogle Scholar
 Helikar T, Rogers JA:ChemChains: a platform for simulation and analysis of biochemical networks aimed to laboratory scientists. BMC Syst Biol. 2009, 3: 58PMID: 19500393 PMCID: PMC2705353,PubMed CentralView ArticlePubMedGoogle Scholar
 Krumsiek J, Pösterl S, Wittmann DM, Theis FJ:Odefy–from discrete to continuous models. BMC Bioinformatics. 2010, 11: 233PMID: 20459647 PMCID: PMC2873544,PubMed CentralView ArticlePubMedGoogle Scholar
 Ferlay J, Shin HR, Bray F, Forman D, Mathers C, Parkin DM:Estimates of worldwide burden of cancer in 2008: Globocan 2008. Int J Cancer. 2010, 127 (12): 28932917.View ArticlePubMedGoogle Scholar
 Henjes F, Bender C, Heyde SVD, Braun L, Mannsperger HA, Schmidt C, Wiemann S, Hasmann M, Aulmann S, Beissbarth T, Korf U:Strong EGFR signaling in cell line models of ERBB2amplified breast cancer attenuates response towards ERBB2targeting drugs. Oncogenesis. 2012, 1 (7): 16View ArticleGoogle Scholar
 Olayioye MA, Neve RM, Lane HA, Hynes NE:The ErbB signaling network: receptor heterodimerization in development and cancer. EMBO J. 2000, 19 (13): PMID: 10880430View ArticleGoogle Scholar
 Heil J, Gondos A, Rauch G, Marmé F, Rom J, Golatta M, Junkermann H, Sinn P, Aulmann S, Debus J, Hof H, Schütz F, Brenner H, Sohn C, Schneeweiss A:Outcome analysis of patients with primary breast cancer initially treated at a certified academic breast unit. Breast (Edinburgh, Scotland). 2012, 21 (3): 303308. PMID: 22310244,View ArticleGoogle Scholar
 Jelovac D, Wolff AC:The adjuvant treatment of HER2positive breast cancer. Curr Treat Options Oncol. 2012, 13 (2): 230239. PMID: 22410709,View ArticlePubMedGoogle Scholar
 Park JW, Neve RM, Szollosi J, Benz CC:Unraveling the biologic and clinical complexities of HER2. Clin Breast Cancer. 2008, 8 (5): 392401. PMID: 18952552,View ArticlePubMedGoogle Scholar
 Tinoco G, Warsch S, Glück S, Avancha K, Montero AJ:Treating breast cancer in the 21st century: emerging biological therapies. J Cancer. 2013, 4 (2): 117132. PMID: 23386910,PubMed CentralView ArticlePubMedGoogle Scholar
 Heyde Svd, Beissbarth T:A new analysis approach of epidermal growth factor receptor pathway activation patterns provides insights into cetuximab resistance mechanisms in head and neck cancer. BMC Medicine. 2012, 10 (1): 43PMID: 22548923,PubMed CentralView ArticlePubMedGoogle Scholar
 Motoyama AB, Hynes NE, Lane HA:The efficacy of ErbB receptortargeted anticancer therapeutics is influenced by the availability of epidermal growth factorrelated peptides. Cancer Res. 2002, 62 (11): 31513158.PubMedGoogle Scholar
 Wilson TR, Fridlyand J, Yan Y, Penuel E, Burton L, Chan E, Peng J, Lin E, Wang Y, Sosman J, Ribas A, Li J, Moffat J, Sutherlin DP, Koeppen H, Merchant M, Neve R, Settleman J:Widespread potential for growthfactordriven resistance to anticancer kinase inhibitors. Nature. 2012, 487 (7408): 505509. PMID: 22763448,PubMed CentralView ArticlePubMedGoogle Scholar
 Gallardo A, Lerma E, Escuin D, Tibau A, Muñoz J, Ojeda B, Barnadas A, Adrover E, SánchezTejada L, Giner D, OrtizMartínez F, Peiró G:Increased signalling of EGFR and IGF1R, and deregulation of PTEN/PI3K/Akt pathway are related with trastuzumab resistance in HER2 breast carcinomas. Br J Cancer. 2012, 106 (8): 13671373. PMID: 22454081,PubMed CentralView ArticlePubMedGoogle Scholar
 Wang L, Zhang Q, Zhang J, Sun S, Guo H, Jia Z, Wang B, Shao Z, Wang Z, Hu X:PI3K pathway activation results in low efficacy of both trastuzumab and lapatinib. BMC Cancer. 2011, 11: 248PMID: 21676217,PubMed CentralView ArticlePubMedGoogle Scholar
 Diermeier S, HorvÃąth G, KnuechelClarke R, Hofstaedter F, Söllosi J, Brockhoff G:Epidermal growth factor receptor coexpression modulates susceptibility to herceptin in HER2/neu overexpressing breast cancer cells via specific erbBreceptor interaction and activation. Exp Cell Res. 2005, 304 (2): 604619. PMID: 15748904,View ArticlePubMedGoogle Scholar
 Pallis AG, Syrigos KN:Epidermal growth factor receptor tyrosine kinase inhibitors in the treatment of NSCLC. Lung cancer (Amsterdam, Netherlands). (2013), PMID: 23384674,Google Scholar
 Moore MJ, Goldstein D, Hamm J, Figer A, Hecht JR, Gallinger S, Au HJ, Murawa P, Walde D, Wolff RA, Campos D, Lim R, Ding K, Clark G, VoskoglouNomikos T, Ptasynski M, Parulekar W, National Cancer, Institute of Canada Clinical Trials Group:Erlotinib plus gemcitabine compared with gemcitabine alone in patients with advanced pancreatic cancer: a phase III trial of the national cancer institute of canada clinical trials group. J Clin Oncol. 2007, 25 (15): 19601966. PMID: 17452677,View ArticlePubMedGoogle Scholar
 Loebke C, Sueltmann H, Schmidt C, Henjes F, Wiemann S, Poustka A, Korf U:Infraredbased protein detection arrays for quantitative proteomics. PROTEOMICS. 2007, 7 (4): 558564.View ArticlePubMedGoogle Scholar
 Kataoka Y, Mukohara T, Shimada H, Saijo N, Hirai M, Minami H:Association between gainoffunction mutations in PIK3CA and resistance to HER2targeted agents in HER2amplified breast cancer cell lines. Ann Oncol. 2010, 21 (2): 255262. PMID: 19633047,View ArticlePubMedGoogle Scholar
 Nahta R, Yuan LXH, Zhang B, Kobayashi R, Esteva FJ:Insulinlike growth factori receptor/human epidermal growth factor receptor 2 heterodimerization contributes to trastuzumab resistance of breast cancer cells. Cancer Res. 2005, 65 (23): 1111811128. PMID: 16322262,View ArticlePubMedGoogle Scholar
 Hornbeck PV, Kornhauser JM, Tkachev S, Zhang B, Skrzypek E, Murray B, Latham V, Sullivan M:PhosphoSitePlus: a comprehensive resource for investigating the structure and function of experimentally determined posttranslational modifications in man and mouse. Nucleic Acids Res. 2011, 40 (D1): 261270.View ArticleGoogle Scholar
 FrÃűhlich H, Sahin O, Arlt D, Bender C, Beissbarth T:Deterministic effects propagation networks for reconstructing protein signaling networks from multiple interventions. BMC Bioinformatics. 2009, 10: 322PMID: 19814779,View ArticleGoogle Scholar
 Bender C, Henjes F, Fröhlich H, Wiemann S, Korf U, Beissbarth T:Dynamic deterministic effects propagation networks: learning signalling pathways from longitudinal protein array data. Bioinformatics (Oxford, England). 2010, 26 (18): 596602. PMID: 20823327,View ArticleGoogle Scholar
 Almudevar A, McCall MN, McMurray H, Land H:Fitting boolean networks from steady state perturbation data. Stat Appl Genet Mol Biol. 2011, 10 (1): 140.Google Scholar
 Brooks SP, Gelman A:General methods for monitoring convergence of iterative simulations. J Comput Graph Stat. 1998, 7 (4): 434455.Google Scholar
 Benjamini Y, Hochberg Y:Controlling the false discovery rate: A practical and powerful approach to multiple testing. J R Stat Soc Series B (Methodological). 1995, 57 (1): 289300. ArticleType: researcharticle/Full publication date: 1995/Copyright Ⓒ1995 Royal Statistical Society,Google Scholar
 Odenbrett MR, Wijs A, Ligtenberg W, Hilbers P, Bo Na Ki, D:Efficient reconstruction of biological networks via transitive reduction on general purpose graphics processors. BMC Bioinformatics. 2012, 13 (1): PMID: 23110660Google Scholar
 Mikalsen T, Gerits N, Moens U:Inhibitors of signal transduction protein kinases as targets for cancer therapy. M Raafat ElGewely (ed.) Biotechnology Annual Review. Elsevier, 153223. ISBN: 13872656 2006. http://www.sciencedirect.com/science/article/pii/S1387265606120062.,
 Dienstmann R, De Dosso S, Felip E, Tabernero J:Drug development to overcome resistance to EGFR inhibitors in lung and colorectal cancer. Mol Oncol. 2012, 6 (1): 1526.View ArticlePubMedGoogle Scholar
 Esteva FJ, Pusztai L:Optimizing outcomes in HER2positive breast cancer: the molecular rationale. Oncology (Williston Park, N.Y.). 2005, 19 (13 Suppl 5): 516. PMID: 19364051,Google Scholar
 Sato S, Fujita N, Tsuruo T:Involvement of 3phosphoinositidedependent protein kinase1 in the MEK/MAPK signal transduction pathway. J Biol Chem. 2004, 279 (32): 3375933767. PMID: 15175348,View ArticlePubMedGoogle Scholar
 Maurer M, Su T, Saal LH, Koujak S, Hopkins BD, Barkley CR, Wu J, Nandula S, Dutta B, Xie Y, Chin YR, Kim DI, Ferris JS, GruvbergerSaal SK, Laakso M, Wang X, Memeo L, Rojtman A, Matos T, Yu JS, CordonCardo C, Isola J, Terry MB, Toker A, Mills GB, Zhao JJ, Murty VVVS, Hibshoosh H, Parsons R:3phosphoinositidedependent kinase 1 potentiates upstream lesions on the phosphatidylinositol 3kinase pathway in breast carcinoma. Cancer Res. 2009, 69 (15): 62996306. PMID: 19602588,PubMed CentralView ArticlePubMedGoogle Scholar
 Tseng PH, Wang YC, Weng SC, Weng JR, Chen CS, Brueggemeier RW, Shapiro CL, Chen CY, Dunn SE, Pollak M, Chen CS:Overcoming trastuzumab resistance in HER2overexpressing breast cancer cells by using a novel celecoxibderived phosphoinositidedependent kinase1 inhibitor. Mol Pharmacol. 2006, 70 (5): PMID: 16887935View ArticleGoogle Scholar
 Vega F, Medeiros LJ, Leventaki V, Atwell C, ChoVega JH, Tian L, Claret FX, Rassidakis GZ:Activation of mammalian target of rapamycin signaling pathway contributes to tumor cell survival in anaplastic lymphoma kinasepositive anaplastic large cell lymphoma. Cancer Res. 2006, 66 (13): 65896597. PMID: 16818631,View ArticlePubMedGoogle Scholar
 Frödin M, Jensen CJ, Merienne K, Gammeltoft S:A phosphoserineregulated docking site in the protein kinase RSK2 that recruits and activates PDK1. EMBO J. 2000, 19 (12): 29242934. PMID: 10856237,PubMed CentralView ArticlePubMedGoogle Scholar
 Klos KS, Wyszomierski SL, Sun M, Tan M, Zhou X, Li P, Yang W, Yin G, Hittelman WN, Yu D:ErbB2 increases vascular endothelial growth factor protein synthesis via activation of mammalian target of rapamycin/p70S6K leading to increased angiogenesis and spontaneous metastasis of human breast cancer cells. Cancer Res. 2006, 66 (4): 20282037. PMID: 16489002,View ArticlePubMedGoogle Scholar
 Schaefer G, Shao L, Totpal K, Akita RW:Erlotinib directly inhibits HER2 kinase activation and downstream signaling events in intact cells lacking epidermal growth factor receptor expression. Cancer Res. 2007, 67 (3): 12281238. PMID: 17283159,View ArticlePubMedGoogle Scholar
 Zhong Q, Simonis N, Li QR, Charloteaux B, Heuze F, Klitgord N, Tam S, Yu H, Venkatesan K, Mou D, Swearingen V, Yildirim MA, Yan H, Dricot A, Szeto D, Lin C, Hao T, Fan C, Milstein S, Dupuy D, Brasseur R, Hill DE, Cusick ME, Vidal M:Edgetic perturbation models of human inherited disorders. Mol Syst Biol. 2009, 5: 321PMID: 19888216,PubMed CentralView ArticlePubMedGoogle Scholar
 Connor TM, Knudson CM, Korsmeyer SJ, Lowe SW, McCurrach ME:baxdeficiency promotes drug resistance and oncogenic transformation by attenuating p53dependent apoptosis. Proc Nat Acad Sci USA. 1997, 94 (6): 23452349. PMID: 9122197,PubMed CentralView ArticlePubMedGoogle Scholar
 Sherr CJ, McCormick F:The RB and p53 pathways in cancer. Cancer Cell. 2002, 2 (2): 103112. PMID: 12204530,View ArticlePubMedGoogle Scholar
 Sithanandam G, Anderson LM:The ERBB3 receptor in cancer and cancer gene therapy. Cancer Gene Therapy. 2008, 15 (7): 413448. PMID: 18404164,PubMed CentralView ArticlePubMedGoogle Scholar
 Lynch DK, Daly RJ:PKBmediated negative feedback tightly regulates mitogenic signalling via gab2. EMBO J. 2002, 21 (12): 7282. PMID: 11782427,PubMed CentralView ArticlePubMedGoogle Scholar
 Chakrabarty A, Sánchez V, Kuba MG, Rinehart C, Arteaga CL:Feedback upregulation of HER3 (ErbB3) expression and activity attenuates antitumor effect of PI3K inhibitors. Proc Nat Acad Sci USA. 2012, 109 (8): 27182723. PMID: 21368164,PubMed CentralView ArticlePubMedGoogle Scholar
 Abrieu A, Dorée M, Fisher D:The interplay between cyclinbcdc2 kinase (MPF) and MAP kinase during maturation of oocytes. J Cell Sci. 2001, 114 (Pt 2): 257267. PMID: 11148128,PubMedGoogle Scholar
 Jirmanova L, Afanassieff M, GobertGosse S, Markossian S, Savatier P:Differential contributions of ERK and PI3kinase to the regulation of cyclin d1 expression and to the control of the G1/S transition in mouse embryonic stem cells. Oncogene. 2002, 21 (36): 55155528. PMID: 12165850,View ArticlePubMedGoogle Scholar
 Searle JS, Li B, Du W:Targeting rb mutant cancers by inactivating TSC2. Oncotarget. 2010, 1 (3): 228232. PMID: 20706560,PubMed CentralView ArticlePubMedGoogle Scholar
 Liu H, Radisky DC, Nelson CM, Zhang H, Fata JE, Roth RA, Bissell MJ:Mechanism of akt1 inhibition of breast cancer cell invasion reveals a protumorigenic role for TSC2. Proc Nat Acad Sci USA. 2006, 103 (11): 41344139. PMID: 16537497,PubMed CentralView ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. 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.