Integrated experimental and model-based analysis reveals the spatial aspects of EGFR activation dynamics

A Phosphoproteomics Data Resource for Systems-level Modeling of Kinase Signaling Networks

Building mechanistic models of kinase-driven signaling pathways requires quantitative measurements of protein phosphorylation across physiologically relevant conditions, but this is rarely done because of the insensitivity of traditional technologies. By using a multiplexed deep phosphoproteome profiling workflow, we were able to generate a deep phosphoproteomics dataset of the EGFR-MAPK pathway in non-transformed MCF10A cells across physiological ligand concentrations with a time resolution of <12 min and in the presence and absence of multiple kinase inhibitors. An improved phosphosite mapping technique allowed us to reliably identify >46,000 phosphorylation sites on >6600 proteins, of which >4500 sites from 2110 proteins displayed a >2-fold increase in phosphorylation in response to EGF. This data was then placed into a cellular context by linking it to 15 previously published protein databases. We found that our results were consistent with much, but not all previously reported data regarding the activation and negative feedback phosphorylation of core EGFR-ERK pathway proteins. We also found that EGFR signaling is biphasic with substrates downstream of RAS/MAPK activation showing a maximum response at <3ng/ml EGF while direct substrates, such as HGS and STAT5B, showing no saturation. We found that RAS activation is mediated by at least 3 parallel pathways, two of which depend on PTPN11. There appears to be an approximately 4-minute delay in pathway activation at the step between RAS and RAF, but subsequent pathway phosphorylation was extremely rapid. Approximately 80 proteins showed a >2-fold increase in phosphorylation across all experiments and these proteins had a significantly higher median number of phosphorylation sites (~18) relative to total cellular phosphoproteins (~4). Over 60% of EGF-stimulated phosphoproteins were downstream of MAPK and included mediators of cellular processes such as gene transcription, transport, signal transduction and cytoskeletal arrangement. Their phosphorylation was either linear with respect to MAPK activation or biphasic, corresponding to the biphasic signaling seen at the level of the EGFR. This deep, integrated phosphoproteomics data resource should be useful in building mechanistic models of EGFR and MAPK signaling and for understanding how downstream responses are regulated.

Ligand-receptor promiscuity enables cellular addressing

In multicellular organisms, secreted ligands selectively activate, or "address," specific target cell populations to control cell fate decision-making and other processes. Key cell-cell communication pathways use multiple promiscuously interacting ligands and receptors, provoking the question of how addressing specificity can emerge from molecular promiscuity. To investigate this issue, we developed a general mathematical modeling framework based on the bone morphogenetic protein (BMP) pathway architecture. We find that promiscuously interacting ligand-receptor systems allow a small number of ligands, acting in combinations, to address a larger number of individual cell types, defined by their receptor expression profiles. Promiscuous systems outperform seemingly more specific one-to-one signaling architectures in addressing capability. Combinatorial addressing extends to groups of cell types, is robust to receptor expression noise, grows more powerful with increases in the number of receptor variants, and is maximized by specific biochemical parameter relationships. Together, these results identify design principles governing cellular addressing by ligand combinations.

PyDREAM: high-dimensional parameter inference for biological models in python

Summary

Biological models contain many parameters whose values are difficult to measure directly via experimentation and therefore require calibration against experimental data. Markov chain Monte Carlo (MCMC) methods are suitable to estimate multivariate posterior model parameter distributions, but these methods may exhibit slow or premature convergence in high-dimensional search spaces. Here, we present PyDREAM, a Python implementation of the (Multiple-Try) Differential Evolution Adaptive Metropolis [DREAM_(ZS)] algorithm developed by Vrugt and ter Braak (2008) and Laloy and Vrugt (2012). PyDREAM achieves excellent performance for complex, parameter-rich models and takes full advantage of distributed computing resources, facilitating parameter inference and uncertainty estimation of CPU-intensive biological models.

Availability and implementation

PyDREAM is freely available under the GNU GPLv3 license from the Lopez lab GitHub repository at http://github.com/LoLab-VU/PyDREAM.

Supplementary information

Supplementary data are available at Bioinformatics online.

A mechanistic pan-cancer pathway model informed by multi-omics data interprets stochastic cell fate responses to drugs and mitogens

Most cancer cells harbor multiple drivers whose epistasis and interactions with expression context clouds drug and drug combination sensitivity prediction. We constructed a mechanistic computational model that is context-tailored by omics data to capture regulation of stochastic proliferation and death by pan-cancer driver pathways. Simulations and experiments explore how the coordinated dynamics of RAF/MEK/ERK and PI-3K/AKT kinase activities in response to synergistic mitogen or drug combinations control cell fate in a specific cellular context. In this MCF10A cell context, simulations suggest that synergistic ERK and AKT inhibitor-induced death is likely mediated by BIM rather than BAD, which is supported by prior experimental studies. AKT dynamics explain S-phase entry synergy between EGF and insulin, but simulations suggest that stochastic ERK, and not AKT, dynamics seem to drive cell-to-cell proliferation variability, which in simulations is predictable from pre-stimulus fluctuations in C-Raf/B-Raf levels. Simulations suggest MEK alteration negligibly influences transformation, consistent with clinical data. Tailoring the model to an alternate cell expression and mutation context, a glioma cell line, allows prediction of increased sensitivity of cell death to AKT inhibition. Our model mechanistically interprets context-specific landscapes between driver pathways and cell fates, providing a framework for designing more rational cancer combination therapy.

Cancer is a complex and diverse disease. Two people with the same cancer type often respond differently to the same treatment. These differences are primarily driven by the fact that two type-matched tumors can possess distinct sets of mutations and gene expression profiles, provoking differential sensitivity to drugs. Over the past few decades, we have seen a shift away from more broadly cytotoxic drugs to more targeted molecules therapies; but how to match a patient with a specific drug or drug cocktail remains a difficult problem. Here, we build a mechanistic ordinary differential equation model describing the interactions between commonly mutated pan-cancer signaling pathways—receptor tyrosine kinases, Ras/RAF/ERK, PI3K/AKT, mTOR, cell cycle, DNA damage, and apoptosis. We develop methods for how to tailor the model to multi-omics data from a specific biological context, devise a novel stochastic algorithm to induce non-genetic cell-to-cell fluctuations in mRNA and protein quantities over time, and train the model against a wealth of biochemical and cell fate data to gain insight into the systems-level, context-specific control of proliferation and death. One day, we hope models of this kind could be tailored to patient-derived tumor mRNA sequencing data and used to prioritize patient-specific drug regimens.

Quantitative investigation of MDA-MB-231 breast cancer cell motility: dependence on epidermal growth factor concentration and its gradient

Enhanced cell motility is one of the primary features of cancer. Accumulated evidence demonstrates that Epidermal Growth Factor Receptor (EGFR) mediated pathways play an important role in breast cancer cell proliferation and migration. We have quantified the MDA-MB-231 breast cancer cell migration in response to the stimulation of EGFR pathways with their ligand EGF to determine how the cell motility of MDA-MB-231 cells depends on the ligand concentration and gradient. Analysis at the single cell level combined with mathematical modeling and the ability to vary the ligand concentration and gradients locally using microfluidic devices allowed us to separate the unique contributions of ligand concentration and ligand gradient to cell motility. We tracked the motility of 6600 cells individually using time lapse imaging under varying EGF stimulation conditions. Trajectory analysis of the tracked cells using non-linear multivariate regression models showed that: (i) cell migration of MDA-MB-231 breast cancer cells depends on the ligand gradient but not on the ligand concentration. This observation was valid for both the total (direction independent) and directed (along gradient direction) cell velocities. Although the dependence of the directed motility on ligand gradient is to be expected, the dependence of the total velocity solely on ligand gradient was an unexpected novel observation. (ii) Enhancement of the motilities of individual cells in a population upon exposure to the ligand was highly heterogeneous, and only a very small percentage of cells responded strongly to the external stimuli. Separating out the non-responding cells using quantitative analysis of individual cell motilities enabled us to establish that enhanced motility of the responding cells indeed increases monotonically with increasing EGF gradient. (iii) A large proportion of cells in a population were unresponsive to ligand stimulation, and their presence introduced considerable random intrinsic variability to the observations. This indicated that studying cell motilities at the individual cell level is necessary to better capture the biological reality and that population averaging methods should be avoided. Studying motilities at the individual cell level is particularly important to understand the biological processes that are possibly driven by the action of a small portion of cells in a population, such as metastasis. We discuss the implications of our results on the total and chemotactic movement of cancer cells in the tumor microenvironment.

EGFR signaling pathways are wired differently in normal 184A1L5 human mammary epithelial and MDA-MB-231 breast cancer cells

Because of differences in the downstream signaling patterns of its pathways, the role of the human epidermal growth factor family of receptors (HER) in promoting cell growth and survival is cell line and context dependent. Using two model cell lines, we have studied how the regulatory interaction network among the key proteins of HER signaling pathways may be rewired upon normal to cancerous transformation. We in particular investigated how the transcription factor STAT3 and several key kinases' involvement in cancer-related signaling processes differ between normal 184A1L5 human mammary epithelial (HME) and MDA-MB-231 breast cancer epithelial cells. Comparison of the responses in these cells showed that normal-to-cancerous cellular transformation causes a major re-wiring of the growth factor initiated signaling. In particular, we found that: i) regulatory interactions between Erk, p38, JNK and STAT3 are triangulated and tightly coupled in 184A1L5 HME cells, and ii) STAT3 is only weakly associated with the Erk-p38-JNK pathway in MDA-MB-231 cells. Utilizing the concept of pathway substitution, we predicted how the observed differences in the regulatory interactions may affect the proliferation/survival and motility responses of the 184A1L5 and MDA-MB-231 cells when exposed to various inhibitors. We then validated our predictions experimentally to complete the experiment-computation-experiment iteration loop. Validated differences in the regulatory interactions of the 184A1L5 and MDA-MB-231 cells indicated that instead of inhibiting STAT3, which has severe toxic side effects, simultaneous inhibition of JNK together with Erk or p38 could be a more effective strategy to impose cell death selectively to MDA-MB-231 cancer cells while considerably lowering the side effects to normal epithelial cells. Presented analysis establishes a framework with examples that would enable cell signaling researchers to identify the signaling network structures which can be used to predict the phenotypic responses in particular cell lines of interest.

Use of mechanistic models to integrate and analyze multiple proteomic datasets

Proteins in cell signaling networks tend to interact promiscuously through low-affinity interactions. Consequently, evaluating the physiological importance of mapped interactions can be difficult. Attempts to do so have tended to focus on single, measurable physicochemical factors, such as affinity or abundance. For example, interaction importance has been assessed on the basis of the relative affinities of binding partners for a protein of interest, such as a receptor. However, multiple factors can be expected to simultaneously influence the recruitment of proteins to a receptor (and the potential of these proteins to contribute to receptor signaling), including affinity, abundance, and competition, which is a network property. Here, we demonstrate that measurements of protein copy numbers and binding affinities can be integrated within the framework of a mechanistic, computational model that accounts for mass action and competition. We use cell line-specific models to rank the relative importance of protein-protein interactions in the epidermal growth factor receptor (EGFR) signaling network for 11 different cell lines. Each model accounts for experimentally characterized interactions of six autophosphorylation sites in EGFR with proteins containing a Src homology 2 and/or phosphotyrosine-binding domain. We measure importance as the predicted maximal extent of recruitment of a protein to EGFR following ligand-stimulated activation of EGFR signaling. We find that interactions ranked highly by this metric include experimentally detected interactions. Proteins with high importance rank in multiple cell lines include proteins with recognized, well-characterized roles in EGFR signaling, such as GRB2 and SHC1, as well as a protein with a less well-defined role, YES1. Our results reveal potential cell line-specific differences in recruitment.

Quantitative analysis reveals how EGFR activation and downregulation are coupled in normal but not in cancer cells

Ubiquitination of the epidermal growth factor receptor (EGFR) that occurs when Cbl and Grb2 bind to three phosphotyrosine residues (pY1045, pY1068 and pY1086) on the receptor displays a sharp threshold effect as a function of EGF concentration. Here we use a simple modelling approach together with experiments to show that the establishment of the threshold requires both the multiplicity of binding sites and cooperative binding of Cbl and Grb2 to the EGFR. While the threshold is remarkably robust, a more sophisticated model predicted that it could be modulated as a function of EGFR levels on the cell surface. We confirmed experimentally that the system has evolved to perform optimally at physiological levels of EGFR. As a consequence, this system displays an intrinsic weakness that causes--at the supraphysiological levels of receptor and/or ligand associated with cancer--uncoupling of the mechanisms leading to signalling through phosphorylation and attenuation through ubiquitination.

Integrated analysis reveals that STAT3 is central to the crosstalk between HER/ErbB receptor signaling pathways in human mammary epithelial cells

Human epidermal growth factor receptors (HER, also known as ErbB) drive cellular proliferation, pro-survival and stress responses by activating several downstream kinases, in particular ERK, p38 MAPK, JNK (SAPK), the PI3K/AKT, as well as various transcriptional regulators such as STAT3. When co-expressed, the first three members of HER family (HER1-3) can form homo- and hetero-dimers, and there is considerable evidence suggesting that the receptor dimers differentially activate intracellular signaling pathways. To better understand the interactions in this system, we pursued multi-factorial experiments where HER dimerization patterns and signaling pathways were rationally perturbed. We measured the activation of HER1-3 receptors and of the sentinel signaling proteins ERK, AKT, p38 MAPK, JNK, STAT3 as a function of time in a panel of human mammary epithelial (HME) cells expressing different levels of HER1-3 stimulated with various ligand combinations. We hypothesized that the HER dimerization pattern is a better predictor of downstream signaling than the total receptor activation levels. We validated this hypothesis using a combination of model-based analysis to quantify the HER dimerization patterns, and by clustering the activation data in multiple ways to confirm that the HER receptor dimer is a better predictor of the signaling through p38 MAPK, ERK and AKT pathways than the total HER receptor expression and activation levels. We then pursued combinatorial inhibition studies to identify the causal regulatory interactions between sentinel signaling proteins. Quantitative analysis of the collected data using the modular response analysis (MRA) and its Bayesian Variable Selection Algorithm (BVSA) version allowed us to obtain a consensus regulatory interaction model, which revealed that STAT3 occupies a central role in the crosstalk between the studied pathways in HME cells. Results of the BVSA/MRA and cluster analysis were in agreement with each other.

Functional aspects of the EGF-induced MAP kinase cascade: a complex self-organizing system approach

The EGF-induced MAP kinase cascade is one of the most important and best characterized networks in intracellular signalling. It has a vital role in the development and maturation of living organisms. However, when deregulated, it is involved in the onset of a number of diseases. Based on a computational model describing a "surface" and an "internalized" parallel route, we use systems biology techniques to characterize aspects of the network's functional organization. We examine the re-organization of protein groups from low to high external stimulation, define functional groups of proteins within the network, determine the parameter best encoding for input intensity and predict the effect of protein removal to the system's output response. Extensive functional re-organization of proteins is observed in the lower end of stimulus concentrations. As we move to higher concentrations the variability is less pronounced. 6 functional groups have emerged from a consensus clustering approach, reflecting different dynamical aspects of the network. Mutual information investigation revealed that the maximum activation rate of the two output proteins best encodes for stimulus intensity. Removal of each protein of the network resulted in a range of graded effects, from complete silencing to intense activation. Our results provide a new "vista" of the EGF-induced MAP kinase cascade, from the perspective of complex self-organizing systems. Functional grouping of the proteins reveals an organizational scheme contrasting the current understanding of modular topology. The six identified groups may provide the means to experimentally follow the dynamics of this complex network. Also, the vulnerability analysis approach may be used for the development of novel therapeutic targets in the context of personalized medicine.

Reconstruction of biofilm images: combining local and global structural parameters

Digitized images can be used for quantitative comparison of biofilms grown under different conditions. Using biofilm image reconstruction, it was previously found that biofilms with a completely different look can have nearly identical structural parameters and that the most commonly utilized global structural parameters were not sufficient to uniquely define these biofilms. Here, additional local and global parameters are introduced to show that these parameters considerably increase the reliability of the image reconstruction process. Assessment using human evaluators indicated that the correct identification rate of the reconstructed images increased from 50% to 72% with the introduction of the new parameters into the reconstruction procedure. An expanded set of parameters especially improved the identification of biofilm structures with internal orientational features and of structures in which colony sizes and spatial locations varied. Hence, the newly introduced structural parameter sets helped to better classify the biofilms by incorporating finer local structural details into the reconstruction process.

Network Modeling Approach to Predict Myofibroblast Differentiation

Fibrotic disease is a major cause of morbidity and mortality and is characterized by the transition of resident fibroblast cells into active myofibroblasts, identified by their expression of alpha smooth muscle actin. Myofibroblast differentiation is regulated by growth factor signaling and mechanical signals transduced through integrins, which converge at focal adhesion proteins (Src and FAK) and MAPK signaling, but lead to divergent outcomes. While details are known about individual pathways, little is known about their interactions. To this end, an ODE-based model of this cell signaling network was developed in parallel with in vitro experiments to analyze potential mechanisms of crosstalk and regulation of αSMA production. We found that cells lacking Src or FAK produce significantly less or more αSMA than wild type cells, respectively. Transforming growth factor beta 1 and fibroblast growth factor signal through ERK and MAPK p38 with different dynamic profiles to increase or decrease αSMA expression, respectively. Our model effectively recreated αSMA expression levels across a set of 22 experimental conditions and matched some features of transient phosphorylation of ERK and p38. These results support a potential mechanism for regulation of fibroblast differentiation: αSMA production is promoted by active p38 and Src and opposed by ERK.

Model-based analysis of HER activation in cells co-expressing EGFR, HER2 and HER3

The HER/ErbB family of receptor tyrosine kinases drives critical responses in normal physiology and cancer, and the expression levels of the various HER receptors are critical determinants of clinical outcomes. HER activation is driven by the formation of various dimer complexes between members of this receptor family. The HER dimer types can have differential effects on downstream signaling and phenotypic outcomes. We constructed an integrated mathematical model of HER activation, and trafficking to quantitatively link receptor expression levels to dimerization and activation. We parameterized the model with a comprehensive set of HER phosphorylation and abundance data collected in a panel of human mammary epithelial cells expressing varying levels of EGFR/HER1, HER2 and HER3. Although parameter estimation yielded multiple solutions, predictions for dimer phosphorylation were in agreement with each other. We validated the model using experiments where pertuzumab was used to block HER2 dimerization. We used the model to predict HER dimerization and activation patterns in a panel of human mammary epithelial cells lines with known HER expression levels in response to stimulations with ligands EGF and HRG. Simulations over the range of expression levels seen in various cell lines indicate that: i) EGFR phosphorylation is driven by HER1-HER1 and HER1-HER2 dimers, and not HER1-HER3 dimers, ii) HER1-HER2 and HER2-HER3 dimers both contribute significantly to HER2 activation with the EGFR expression level determining the relative importance of these species, and iii) the HER2-HER3 dimer is largely responsible for HER3 activation. The model can be used to predict phosphorylated dimer levels for any given HER expression profile. This information in turn can be used to quantify the potencies of the various HER dimers, and can potentially inform personalized therapeutic approaches.

A family of cell surface molecules called the HER receptor family plays important roles in normal physiology and cancer. This family has four members, HER1-4. These receptors convert signals received from the extracellular environment into cell decisions such as growth and survival – a process termed signal transduction. In particular, HER2 and HER3 are over-expressed in a number of tumors, and their expression levels are associated with abnormal growth and poor clinical prognosis. A key step in HER-mediated signal transduction is the formation of dimer complexes between members of this family. Different dimer types have different potencies for activating normal and aberrant responses. Prediction of the dimerization pattern for a given HER expression level may pave the way for personalized therapeutic approaches targeting specific dimers. Towards this end, we constructed a mathematical model for HER dimerization and activation. We determined unknown model parameters by analyzing HER activation data collected in a panel of human mammary epithelial cells that express different levels of the HER molecules. The model enables us to quantitatively link HER expression levels to receptor dimerization and activation. Further, the model can be used to support additional quantitative investigations into the basic biology of HER-mediated signal transduction.