Cell signaling heterogeneity is modulated by both cell-intrinsic and -extrinsic mechanisms: An integrated approach to understanding targeted therapy

During the last decade, our understanding of cancer cell signaling networks has significantly improved, leading to the development of various targeted therapies that have elicited profound but, unfortunately, short-lived responses. This is, in part, due to the fact that these targeted therapies ignore context and average out heterogeneity. Here, we present a mathematical framework that addresses the impact of signaling heterogeneity on targeted therapy outcomes. We employ a simplified oncogenic rat sarcoma (RAS)-driven mitogen-activated protein kinase (MAPK) and phosphoinositide 3-kinase-protein kinase B (PI3K-AKT) signaling pathway in lung cancer as an experimental model system and develop a network model of the pathway. We measure how inhibition of the pathway modulates protein phosphorylation as well as cell viability under different microenvironmental conditions. Training the model on this data using Monte Carlo simulation results in a suite of in silico cells whose relative protein activities and cell viability match experimental observation. The calibrated model predicts distributional responses to kinase inhibitors and suggests drug resistance mechanisms that can be exploited in drug combination strategies. The suggested combination strategies are validated using in vitro experimental data. The validated in silico cells are further interrogated through an unsupervised clustering analysis and then integrated into a mathematical model of tumor growth in a homogeneous and resource-limited microenvironment. We assess posttreatment heterogeneity and predict vast differences across treatments with similar efficacy, further emphasizing that heterogeneity should modulate treatment strategies. The signaling model is also integrated into a hybrid cellular automata (HCA) model of tumor growth in a spatially heterogeneous microenvironment. As a proof of concept, we simulate tumor responses to targeted therapies in a spatially segregated tissue structure containing tumor and stroma (derived from patient tissue) and predict complex cell signaling responses that suggest a novel combination treatment strategy.

activated protein kinase (MAPK) and phosphoinositide 3-kinase-protein kinase B (PI3K-AKT) signaling pathway in lung cancer as an experimental model system and develop a network model of the pathway. We measure how inhibition of the pathway modulates protein phosphorylation as well as cell viability under different microenvironmental conditions. Training the model on this data using Monte Carlo simulation results in a suite of in silico cells whose relative protein activities and cell viability match experimental observation. The calibrated model predicts distributional responses to kinase inhibitors and suggests drug resistance mechanisms that can be exploited in drug combination strategies. The suggested combination strategies are validated using in vitro experimental data. The validated in silico cells are further interrogated through an unsupervised clustering analysis and then integrated into a mathematical model of tumor growth in a homogeneous and resource-limited microenvironment. We assess posttreatment heterogeneity and predict vast differences across treatments with similar efficacy, further emphasizing that heterogeneity should modulate treatment strategies. The signaling model is also integrated into a hybrid cellular automata (HCA) model of tumor growth in a spatially heterogeneous microenvironment. As a proof of concept, we simulate tumor responses to targeted therapies in a spatially segregated tissue structure containing tumor and stroma (derived from patient tissue) and predict complex cell signaling responses that suggest a novel combination treatment strategy. PLOS

Introduction
Normal cell signaling is significantly altered in cancer as a result of genetic and epigenetic changes, facilitating uncontrolled proliferation and cell survival [1,2]. Targeted therapies directly exploit these alterations by blocking the activity of specific proteins typically mutated or abnormally up-regulated [3]. These therapies have elicited dramatic success in controlling the growth of multiple cancers [4][5][6][7][8][9][10] but showed little to moderate impact on others [11,12]. Drug resistance, however, remains a major problem due to both cancer cell-intrinsic (innate and acquired) resistance mechanisms [13] and microenvironment-mediated resistance [14][15][16]. Tumor heterogeneity is known to contribute to drug resistance [17,18]. Cancer cells within a tumor exhibit differential genetic and phenotypic characteristics [19]. Genomic heterogeneity leads to cell-to-cell variability in protein expression and activity as genes drive the production of proteins. Protein activity is variable even in genetically identical cancer cell populations in the same microenvironment [20][21][22]. This cell-to-cell variability arises from intrinsic stochastic fluctuations [23][24][25][26][27][28][29][30] and variation in microenvironmental conditions that affect the protein-signaling network. This variation can affect sensitivity to stimuli, contribute to cell phenotype decisions, and cause clonal cells to differently respond to stimulus and targeted therapies (e.g., erlotinib) [31,32].
Understanding how cancer cell signaling variation affects targeted therapy outcomes is challenging. Cancer-driven signaling proteins do not function in isolation but rather function in protein complexes that belong to large and complex signaling networks that govern key phenotypic processes such as proliferation, apoptosis, and response to targeted therapy [33,34]. Furthermore, the cancer signaling network and drug response are modulated by microenvironmental factors [35,36]. Therefore, experimental data obtained from simplified cell-based experiments in single uniform environments will have limited ability to tease apart the impact of signaling network and microenvironmental variation on targeted therapy outcomes. A number of previous studies have used mathematical models to understand complex signaling networks and improve treatment strategies. Various modeling approaches were developed (reviewed in [37]). Boolean or probabilistic Boolean models were developed to analyze cancer signal pathways and predict treatment outcomes [38][39][40][41]. A logical modeling approach was employed to understand various cell signaling pathways [42][43][44][45][46][47]. An artificial neural network approach was used to map between microenvironments, pathways, and phenotypes [48]. Detailed kinetic models of cell signaling pathways have been studied using systems of ordinary differential equations (ODEs) [49][50][51][52][53][54][55][56][57]. Most studies, however, ignore signaling heterogeneity or extrinsic variation in microenvironmental cues that will differentially stimulate the signaling network.
To investigate the effects of signaling heterogeneity on targeted therapy outcomes, we develop an integrated approach combining in vitro experiments with three different mathematical models, an intracellular signaling model, a cancer cell population growth model, and a hybrid cellular automata (HCA) model of tumor and stroma (see Fig 1 for an overview).

Mathematical modeling of MAPK and PI3K-AKT pathway
We consider a simplified signaling network composed of interactions between key proteins in an oncogenic rat sarcoma (RAS)-driven mitogen-activated protein kinase (MAPK) and phosphoinositide 3-kinase (PI3K)-protein kinase B (AKT, PKB) pathway (Fig 2A). This pathway has been studied extensively, and the positive or negative feedback regulations between proteins in the pathway are known [58][59][60][61][62]. Of note, mammalian cells express three RAS gene family members (HRAS, NRAS, and KRAS), and our model is based on empirical data obtained from a lung cancer cell line (A549); KRAS is always activated by a point mutation (Gly12Ser), whereas the other two RAS proteins (HRAS and NRAS) are wild type. Recent studies reported the importance of crosstalk between wild-type and mutant RAS proteins in cancers driven by oncogenic mutant RAS [59,63,64]. Therefore, we consider two different types of RAS, wild type (RAS_w) and mutant type (RAS_m). The network node connectivity is based on prior pathway information between signaling proteins [58][59][60][61][62]. Most interactions are feed-forward and positive (Fig 2A, green lines) except the one negative feedback regulation of epidermal growth factor receptor (EGFR) by extracellular receptor kinase (ERK) (Fig 2A, red  line). While network connectivity is assumed fixed in the model, the strength of interactions is variable and is modeled using a weight matrix (W). Each element in the network (node x i ) is updated by solving the following equation, . . . ; n; where x 1,2 correspond to two inputs (growth factor and hepatocyte growth factor [HGF]), and χ 3,4. . .,n correspond to the relative change of protein activities or cell viability due to an inhibitor with respect to untreated conditions (log 2 (treated/untreated)) to be consistent with experimental data. Of note, all of the experimental measures in our study are relative values, normalized to the unperturbed condition. The absolute concentration or activity of signaling proteins as well as cell viability are difficult to acquire from experiments performed in our study, namely western blot experiments and cell viability measurement assays. Therefore, the weights in the model represent relative abundance or protein activities in treated conditions compared to treatment-naïve conditions. The model assumes that the rate of change of a variable is determined by the linear combination of neighboring nodes with corresponding weights. This additive linear function has successfully described protein reaction networks [54,55,65] although other functions such as Michaelis-Menten kinetics are viable options [51]. In the experiments we carried out, the microenvironmental conditions are growth factor and HGF. The growth factor (model variable x 1 ) is always present, while HGF (model variable x 2 ) is present in only some of the experimental conditions. In particular, the HGF is not present in the control condition. All experimental results are normalized to this control (no-HGF) condition. To represent these experimental conditions in the model, the input value x 1 is set to a nonzero constant value (e.g., x 1 C), while the variable x 2 is set to be 0 or nonzero (e.g., x 2 C) if HGF is present. We chose to set the parameter constant C to be 10. N(i) represents the neighborhood of a protein node i (a set of nodes connected to the node i), and α indicates a tendency to return to the untreated state. The transfer function T accounts for saturation effects, and the constants ε and β modulate amplitude and slope. In the model, we set ε to be 4.5 and β to be 0.5 to model a smooth sigmoidal behavior. Drug inhibition is modeled by knocking out a corresponding protein activity. For example, if a drug inhibits a protein x i , then the variable x i is set to be a very small number (x i ¼ log 2 x treated x untreated ¼ À m; m ) 1). It is worth noting that we modeled the effect of drugs tyrosineprotein kinase Met inhibitor (METi), EGFRi, and AKTi by inhibiting pMET, pEGFR, and pAKT activities, respectively. Two other drugs, mitogen-activated protein kinase kinase (MEKi), and ERKi, are modeled to inhibit pERK and pRSK, respectively ( Fig 2B).

Logistic tumor growth
We combine the pathway model to a logistic growth ODE model to simulate the growth of in silico cells as a well-mixed population in a resource-limited environment. The model is defined as follows, where y i represents the number of in silico cell i, r i is the cell intrinsic growth rate of the cell i, N is the total number of cell types, and K is the carrying capacity (set to be 1 billion). To model influence of the signaling pathway on cell population growth, we formulate a cell population  . Specifically, the cell viability for each cell type is obtained by fitting pathway model to our experimental data, as described in the Results section. The cell viability of cell type i (θ i ) represents the number of cell type i that survived after being given therapy relative to an untreated control condition (i.e., θ i = P i (t)/P 0 (t), θ i : cell viability, P i,0 (t): number of cell type i at time t in a treated and untreated condition, respectively). To obtain a functional form for the growth rate, we make the following assumptions. We assume that a cell population initially grows exponentially ðP i ðtÞ ¼ P i ð0Þe r i t ; P 0 ðtÞ ¼ P 0 ð0Þe r 0 t ; t < T for some time T). We also assume that the number of initial cell population is the same in a treated condition and in an untreated control condition (P i (0) = P 0 (0)). Then, cell viability is presented as a function of growth rates and time (y i ¼ P i ðtÞ=P 0 ðtÞ ¼ e r i t =e r 0 t ). Solving the function for growth rate r i , we obtain a functional form for the growth rate. We use the doubling time of A549 cells (22 hours from [66]) to obtain the treatment-naïve growth rate (r 0 = 0.76 per day). We use our cell viability assay experimental time point (t = 3 days, described in Results section). Now, we have a constant growth rate of cell type i for each treatment condition (500 cells x 28 treatment conditions, total 14,000 growth rates, r i ). Then, the ODE Eq (2) is solved to simulate a given treatment response of cell i over time. All of the in silico cells are solved simultaneously competing for limited resource (carrying capacity K).

HCA tumor model
The pathway model is integrated into an HCA model [67,68] to simulate treatment responses in a spatially heterogeneous microenvironment. The model has the following assumptions. In the HCA, the cells are defined as points on a two-dimensional lattice that also contains continuous concentration fields of microenvironmental factors, together representing a cross-section of tumor composed of cancer cells and stroma (50 cells x 38 cells). Here, we define the tumor and stroma region explicitly based on an image segmentation of lung adenocarcinoma tissue from a patient. The tumor region contains cancer cells. Each cancer cell contains the pathway model, as developed above (Fig 2A), that links the microenvironment to cell phenotypes. The model grid can contain any number of possible microenvironmental variables. For simplicity, however, we consider only growth factors and HGF. The growth factors are assumed to be constant in the domain. We explicitly model HGF dynamics in space and time using the following partial differential equation, @Hðx; tÞ @t ¼ Dr 2 Hðx; tÞ À lHðx; tÞ; Erlotinib), RAF inhibitor (RAFi, LY3009120), MEK inhibitor (MEKi, GDC0623), ERK inhibitor (ERKi, SCH772984), and AKT inhibitor (AKTi, MK2208) in both control medium (DMSO) and after 2-hour stimulation-by-HGF (50 ng/mL) condition (DMSO plus HGF). (C) Relative cell viabilities after treatments. Cells were treated inhibitors (1 μM) for 72 hours. Cell viabilities were assessed by CellTiter-Glo assay (Promega). Representative triplicates (± SD) are presented, which showed similar results at least three times. (D) The western blots were quantified using ImageJ and relative changes (log2 scale) are reported. Average values of relative cell viabilities are also reported in log2 scale. All the data are normalized to the treatment-naïve control condition. pAKT (T308) readouts are quantified and used in the model. Of note, we didn't quantify total protein levels because our primary interest was protein activity (protein phosphorylation where H(x,t) represents the concentration of HGF at a lattice point x in tumor region and at time t. D represents the diffusion rate, and λ is a decay rate. The parameter values used in a simulation are given in the corresponding figure legend. The concentration of HGF is fixed to be a constant value (H(x,t) = γ) in the stromal region. A Neumann boundary condition ( @H @n x ð Þ ¼ 0, normal derivative = 0) was imposed on the domain boundary. A Dirichlet condition (H(x,t) = γ) was imposed at the interface between tumor and stroma.
The steady state solution of Eq (3) is fed into one of the inputs (HGF) in the pathway model in each cell (Fig 2A). The pathway determines cell viability and controls three different phenotypes-proliferation, quiescence, and death-as defined by the rules summarized in the flowchart (S1 Fig). Each cell is allowed to execute only one phenotype per time step (day). Of note, the model considers only orthogonal neighbors (north, west, south, and east) for space to divide or move. Cells are not allowed to leave nor enter across the boundary and are thus confined within the domain. We assume that the distribution of HGF barely changes during HCA simulation. For example, if a cell divides into two daughter cells, this increased number of cells does not impact the HGF distribution because we assume the consumption of HGF by cells negligible and that cells do not produce HGF. Because HGF consumption is minimal and the HGF diffusion timescale (approximate seconds) is a lot shorter than cell division time scale (approximate day) and the model domain size is small (50 cells x 38 cells), any consumption would quickly be equilibrated to the steady state.

Cell line
A549 lung adenocarcinoma cell line was maintained in RPMI 1640 medium supplemented with 10% FBS. Cells were confirmed to be free of mycoplasma using PlasmoTest (Invivogen, San Diego, CA).

Western blots
Cells were washed with ice-cold PBS, and whole cell extracts were prepared using lysis buffer (0.5% NP-40, 50 mM Tris-Cl, pH 8.0, 150 mM NaCl, 1 mM EDTA) supplemented with protease inhibitor (Roche, Mannheim, Germany) and phosphatase inhibitor cocktail (Sigma-Aldrich, Carlsbad, CA). Whole-cell extracts were resolved on SDS-PAGE and transferred to nitrocellulose membrane. The membrane was blocked in 5% skim milk/PBST and then incubated in primary antibody at 4˚C overnight. Bound antibodies were visualized by horseradish peroxidase-conjugated secondary antibodies and SuperSignal West Pico Chemiluminescent Substrate (Thermo Scientific, Waltham, MA). Primary antibodies used for our study were purchased from Cell Signaling Technology (Danvers, MA) (except for β-actin, which was from Sigma-Aldrich, St. Louis, MO).

Cell viability measurement
Cells were plated on 96-well plate at 2,000 cells per well and then exposed to drugs for 72 hours. Cell viability was analyzed by CellTiter-Glo (Promega, Madison, WI) according to the manufacturer's recommendations.

Proximity ligation assays
Proximity ligation assays were performed as described using Duolink Far Red kit (Sigma-Aldrich, Carlsbad, CA) with antibodies to the following: EGFR (clone B38, Cell Signaling, Danvers, MA), GRB2 (clone 81, BD, San Jose, CA), and AF488-conjugated pan cytokeratin (clone Ae1/Ae3, eBioSciences, San Diego, CA). Tissue specimen was from a de-identified patient treated at Moffitt Cancer Center and was obtained via an institutionally approved protocol.

Results
We followed the steps described in Fig 1 to investigate the effects of signaling heterogeneity on targeted therapy outcomes. First, we develop an intracellular signaling pathway model. We construct a simplified cancer signaling pathway based on prior information about the pathway, and we experimentally perturb the pathway using various kinase inhibitors in two different microenvironmental conditions (Fig 1 Design & Experiments). Among key microenvironmental factors, HGF has been shown to contribute to resistance in multiple Food and Drug Administration (FDA)-approved targeted therapy drugs [35,36]. Therefore, we consider HGF as an additional microenvironmental stimulus. Then, we build a mathematical model of the cancer signaling pathway to predict both signaling and a phenotypic response (cell viability change due to a given therapy) to different inhibitors that target the pathway (Fig 1 prior pathway information). It is important to note that the model includes two microenvironmental factors (i.e., growth factor or HGF) and cell viability ( contains the trained network model that links microenvironment to its phenotype, which determines cell fate in a given condition. As a proof of concept, we simulate the response to an inhibitor in a section of tissue composed of both tumor cells and stroma. The model predicts complex cancer cell signaling responses and treatment outcomes, driven by both cell-intrinsic and -extrinsic mechanisms.

Experimentally measured cancer cell response to various kinase inhibitors
Kirsten rat sarcoma (KRAS)-driven cancer treatment is an important clinical need that remains largely unmet due to limited targeted drug efficacy of key downstream effectors, including MAPK and PI3K-AKT pathways. We therefore choose the KRAS mutant non-small-cell lung cancer (NSCLC) cell line (A549 cell line) as our experimental model system. Using our simplified oncogenic KRAS signaling pathway (Fig 2A, Mathematical modeling section for an explanation of how we obtained this simplified pathway) as a guide, we pharmacologically inhibited individual proteins in the MAPK pathway (MET, EGFR, MEK, ERK inhibitors) and PI3K-AKT (AKT inhibitor) pathway in A549 cells in the both absence and presence of HGF. Drug-induced changes in the phosphorylation of pathway proteins, surrogates for protein activity, were measured by western blotting (Fig 2B). Cell viability was also assessed after 72 hours of drug treatment (Fig 2C). These experimental data were quantified using ImageJ ( Fig  2D), and all data were normalized to the control experimental condition (treatment-naïve condition).

Model calibration
The quantified changes ( Fig 2D) and the Monte Carlo simulation were then employed as an optimization procedure to estimate model parameters (weights, W ij ) that minimize our cost function. A network with lower cost represents the experimental data more accurately. The cost function is defined as follows, where x i;d is the steady state activity of protein or cell viability x i in treatment condition d, y d represents experimental data, and M is the total number of treatment conditions. The weight w jr indicates the weight between RAS_m (r) and its neighbors (N(r)), and the constants (χ, η) modulate the magnitude of the penalty. The first term explains the difference between model prediction and experimental data for a network W. The second term is directed at the activating RAS mutation and incorporates a penalty for estimated weights from the RAS_m node (mutant RAS) that are too small. We included this penalty because our model is based on empirical data of a KRAS mutant cancer cell line (A549 cell), where the resulting KRAS protein is constitutively active. We aimed to capture this activating mutation by penalizing small weights from RAS_m to its neighbors. We used the following method to implement Monte Carlo simulations: 1. Initialize a sparse weight matrix (W, W ij = 0, for no connection in Fig 2A) with random numbers.
2. Enforce the weight elements to satisfy the prior pathway information (W ij = |W ij | for green line; W ij = −|W ij | for red line; Fig 2A).
7. If C new < C old , accept the perturbation; otherwise, mostly reject the perturbation. Only accept the perturbation with a small probability of p, where p = exp(−ν(|C new − C old |) for some large ν.
8. Go to the step 5 and repeat 5 through 7 until achieving error = |C new − C old | < δ, where δ is a predefined tolerance (for small non-negative number, δ > 0).
The model calibration resulted in more than 5,000 weight matrices that fit to the experimental data. We selected the best 500 (top 10%) weight matrices and used these to define our 500 in silico cells. The distributions of in silico cells are presented as box plots in Fig 2E along with the experimental measures. Errors (root-mean-squared-error [RMSE] formula given in S1 Text) are in the range of (0.03-0.56, except ERK: 0.96). The fit of ERK was poor because of unexpected inhibition of pERK by the drug (ERK inhibitor, SCH772984) [69]. The trained networks (weights) are quite heterogeneous (S1 Table). The distribution for each weight is different (S1 Table, skewed, normal, bimodal distributions, with a range of heterogeneity [Shannon] index values).
The weights here may represent relative protein abundance or protein-binding activity. There is ample evidence for differential abundance of protein species across cellular populations. An excellent example was recently published showing that variations in adaptor protein abundance are a major source of regulation of the EGFR-MAPK pathway [70]. There are several examples of differential binding activity of proteins in cell signal transduction. It is well established that adaptors such as GRB2, SHC1, and GAB1 can be recruited to receptor tyrosine kinases (RTKs) either directly or indirectly. Therefore, stochastic variation in multiprotein complex composition at individual receptors exists, and this will vary both within and between cells. It is also accepted that activation-induced receptor degradation and phosphatase activity will affect not only RTK adaptor interactions but also downstream signaling molecules such as rapidly accelerated fibrosarcoma (RAF), MEK, and ERK. Additionally, we have previously shown that, in EGFR mutant cell lines, only a fraction of the receptor is phosphorylated, and cell lines harboring the same oncogenic mutation have different levels of phosphorylated Tyrosine (Tyr) and Serine/threonine (Ser/Thr) residues [71]. Collectively, these examples indicate that protein-protein interactions in response to growth factors are not simply on-off states and that multiple factors independent of protein abundance control final signaling output.

Distributional responses to kinase inhibitors
We simulated responses of in silico cells to seven different inhibitors (EGFRi, METi, RAS_mi, AKTi, RAFi, MEKi, and ERKi). Of note, cell viability in untreated conditions is set to be a single constant value (cell viability untreated 1). It is also worth noting that RAS_mi is assumed to inhibit only RAS_m (RAS mutant), not RAS_w. Distributions of relative cell viability (log2 scale, log 2 (treated/untreated)) of all in silico cells are presented in Fig 3A. Similar to the experimental results (Fig 2B-2D), MEKi and ERKi reduced average cell viability significantly, whereas mean effects of EGFRi, METi, and AKTi are marginal. These results reveal quite heterogeneous responses to drugs, which could be assessed by experimental approaches [72,73]. For example, the distribution of EGFRi treatment is bimodal due to a bimodal distribution of trained MEK-ERK weights (S1 Table and

In silico-predicted rational drug combinations
With our calibrated in silico cell lines, we next examine which drug combinations significantly reduce cell viability and which show marginal effects. For example, what should be cotargeted with AKTi to decrease cell viability significantly? The model predicts that activation of alternative pathways under therapy may provide an escape route to therapy. For example, AKT inhibitor induced increased activity of ERK and ribosomal S6 kinase (RSK) (S3 Fig). Cells with high ERK and RSK activity display resistance (cell viability >0) under the inhibitor, which implies that simultaneous inhibition of this alternative pathway would overcome resistance. We reasoned that combination targeting of proteins that are highly correlated with relative cell viability under a given treatment would be beneficial. In order to identify such protein nodes, scatter plots between predicted protein activity (phosphorylation) and predicted relative cell viability were considered (S4 Fig). Then, the Pearson's coefficient was calculated for all cases. We observed that, under multiple treatment conditions-including inhibition of EGFR, MET, RAS_m, AKT, and RAF-ERK and RSK showed the highest correlation with relative cell viability (S4 Fig, pink box). This suggests that co-inhibition of ERK (by MEKi) or RSK (by ERKi) activity with other therapies would decrease cell viability most significantly. Further simulation revealed an additional application of either MEKi or ERKi to each-EGFRi, METi, RAS_mi, RAFi, and AKTi-significantly decreased cell viability compared to monotherapy of EGFRi, METi, RAS_mi, RAFi, and AKTi (Fig 3B). We tested some of these model predictions experimentally and validated, to some degree, the model's predictive ability (Fig 3C).

Hierarchical clustering of in silico cell responses
We next systematically assessed cell viability reduction to all mono and combination therapies using an unsupervised hierarchical clustering approach to classify cell populations on the basis of their treatment response (relative cell viability). The treatments were categorized into multiple groups (Fig 4A, a tree diagram on the right end of heat-map). The combination of AKTi with MEKi is uniformly effective to all the cells (see the red asterisk [ Ã ] row, dark blue across all in silico cells, with little variation between cells). This combination (MEKi/AKTi) has previously been shown to be effective in NSCLC both in vitro and in vivo [74]. A striking variation is observed in response to the treatments of AKTi, RAS_mi, and RAS_mi/AKTi (Fig 4A, red bars in the first two groups versus dark blue to yellow bars in the rest of the groups). The first two clusters (pink and black color on the top of the heat-map) are associated with poor responses (little to no reduction of cell viability after a given therapy), while others are correlated with good treatment outcomes (significant reduction of cell viability after a given therapy).

Possible mechanisms of response and resistance
Why are some cells sensitive to a given therapy while others are resistant to the same therapy? We hypothesize that this differential drug sensitivity, at least within the context of our model, must be attributed to differential protein activity as modulated by protein-protein interactions (i.e., the weights). To highlight possible mechanisms, we visualized the weights between protein nodes (W ij in our model) using both circular chord diagrams [75] and network diagrams with weighted edges (Fig 4, bottom panels; left: chord diagram; right: weighted network with different edge widths). In the circular diagrams, protein nodes are arranged around a circle with the weight between protein nodes connected to each other through the use of arcs. The width of each arc is determined proportionally by the weight between two protein nodes. To illustrate differences in signaling, we selected two representative in silico cells (Fig 4A and 4B), where cell a is resistant to AKTi, RAS_mi, and RAS_mi/AKTi and cell b is more sensitive to these therapies. In addition, we compared ranges of all weights of all in silico cells between clusters defined by the hierarchical clustering (S5A Fig). We observe heterogeneity of weights within each cluster and between clusters. Differences between clusters are significant for some weights such as weights of growth factor EGFR, EGFR-RAS_w, MET-RAS_w, RAS_w-RAF, and MEK-ERK (S5A Fig).

Heterogeneous responses to HGF stimulation
We next asked how an additional microenvironmental stimulation would modulate the responses to targeted treatments. To address this question, all in silico cells were treated in the presence of HGF, a significant stromal factor that contributes to drug resistance [35,36]. An unsupervised hierarchical clustering on the basis of cell viability changes, from the no-HGF condition, separated the treatments into several groups (Fig 5A, clustering of treatments on the right tree diagram). The analysis also classified the in silico cells into several groups based on cell viability changes due to HGF stimulation (Fig 5A, tree diagram on the top of heatmap).
To test some of these model predictions, we considered three different combination therapies, AKTi/MEKi, EGFRi/MEKi, and EGFRi/AKTi. Of note, the model predicted that the treatments AKTi/MEKi and EGFRi/MEKi are not affected by HGF stimulation (Fig 5B, first two red bar graphs). The experimental data matched well with these predictions (Fig 5B, first Chord diagram: each node in the circle represents each protein node in the network model, represented by different colors. The thickness of chord between two protein nodes represents a weight between two protein nodes (weight, w ij ). The chords are directed, colored by originating sector color. For example, the interaction between RAS_m and RAF is depicted as a light blue chord because the direction is from RAS_m to RAF (RAS_m ! RAF; color of RAS_m sector: light blue). Weighted network diagram: the width of each edge represents the weight. A thicker edge represents a larger weight between two proteins. The numerical data used in Fig 4 are included in the third sheet S1 Data. AKT (PKB), protein kinase B; EGFR, epidermal growth factor receptor; ERK, extracellular receptor kinase; HGF, hepatocyte growth factor; MEK, mitogen-activated protein kinase kinase; MET (c-MET), tyrosine-protein kinase Met or hepatocyte growth factor receptor (HGFR); PI3K, phosphoinositide 3-kinase; RAF, rapidly accelerated fibrosarcoma; RAS, rat sarcoma; RAS_m, mutated RAS; RAS_w, wild-type RAS; RSK, ribosomal S6 kinase. two gray bar graphs). The model also predicted that the effect of combination therapy EGFRi/ AKTi is significantly modulated by HGF stimulation (Fig 5B, the third red bar). This was also corroborated by experiment ( Fig 5B, the third gray bar).

Possible mechanisms explaining HGF modulation
The suite of in silico cells is differentially affected by HGF stimulation. For example, some cells are not affected by the HGF stimulation (e.g., gray bars in EGFRi/RAS_mi in Fig 5A), while others are significantly affected by stimulation (e.g., red bars in EGFRi/RAS_mi in Fig 5A). Why are some cells affected by HGF stimulation, while others are not? To understand why this is the case, we selected two representative cells (a and b) and visualized the weights between protein nodes using both chord diagrams and weighted network diagrams (Fig 5C, bottom  panels). In cell a, the influence of MET on RAS_w is relatively small (Fig 5C, thin blue chords from MET ! RAS_w). In contrast to cell a, the influence of MET on RAS_w in cell b is much Impact of heterogeneity on targeted therapy stronger ( Fig 5C, cell b, thick blue chord from MET ! RAS_w). This may explain why cell b increased its viability significantly upon HGF stimulation compared to no-HGF condition. In addition, distributions of all weights in each cluster show differential activity of proteins and its effects across all six different clusters (S5B Fig).

Competition between in silico cells drives different dominant populations
So far, treatment responses were assessed as if all in silico cells were treated in isolation. To address the effects of cell competition on treatment outcomes, we simulated all cells growing together under treatment with various inhibitors in a homogeneous, resource-limited microenvironment (Methods section). The entire in silico cell population responds in a similar way on some therapies (Fig 6, e.g., EGFRi, METi, RAFi, EGFRi/RAS_mi, etc.), but under other therapies (Fig 6, e.g., MEKi, ERKi, MEKi/ERKi, AKTi/RAFi, etc.)-due to differential viability -some cells became dominant. Interestingly, some treatments simultaneously selected for the same dominant in silico cell (Fig 6, color-shaded boxes). Addition of HGF to some treatments (e.g., EGFRi, RAS_mi, RAFi) significantly changes the fitness of in silico cells and thus drives a different dominant in silico cell (S6 Fig, dotted lines in no-HGF vs solid lines in HGF).

Assessment of post-treatment in silico cell population heterogeneity
We assessed post-treatment population heterogeneity by measuring the Shannon index (H(x) = −∑ i p i (x)log(p i (x)), where p i (x) is the probability of finding an in silico cell i after a given therapy). We compared the index with average cell viability change (Fig 6B). The two are linearly correlated (ρ = 0.65), implying that the less effective a treatment is in controlling the average in silico population growth, the more heterogeneous the post-treatment population would be (Fig 6B). Combination therapies not only display a better average treatment response but also a less diverse post-treatment population (Fig 6B, boxes versus circles). Among all combination therapies, those that combined either with ERKi or MEKi display much better average therapeutic responses (Fig 6B, small cell viability). These treatments not only effectively decrease average cell viability but also lead to a less diverse post-treatment population (e.g., Fig  6B, EGFRi/[ERKi or MEKi] vs EGFRi/[METi,AKTi,RAFI] in red-color circles). HGF stimulation minimally affected the linear relationship between post-treatment heterogeneity and average cell viability reduction (Fig 6C, ρ = 0.63 vs Fig 6B, ρ = 0.65). However, a few treatments did elicit significant changes in both average response and heterogeneity due to HGF stimulation ( Fig 6D).

Impact of microenvironmental heterogeneity on targeted therapy outcomes
Because activated receptors require multiple protein interactions to activate downstream signaling, we have utilized proximity ligation assays that measure the functional association of RTKs and adaptor proteins. Using NSCLC patient specimens and xenograft models, we have previously identified an association of EGFR:GRB2 complexes and response to EGFR inhibition [76] and more recently identified a correlation between MET:GRB2 complexes and response to MET kinase inhibitors [77]. Using this approach, we consistently observe spatial heterogeneity in abundance of RTKs binding to adaptor complexes. The abundance is often highest at tumor regions that are adjacent to stromal regions (Fig 7A).
Motivated by this experimental observation (Fig 7A), we developed an HCA model to investigate the effects of microenvironmental heterogeneity on treatment outcomes (See Methods section and S1 Fig). A steady state configuration of HGF is considered throughout the whole simulation (See Methods section and Fig 7B). We randomly initialized in silico cells that contain calibrated signaling networks (Fig 2) in the domain to mimic a slice of tumor tissue (Fig 7C first, time step = 0; domain size: 50 cells × 38 cells). As proof of concept, we simulated RAS_m inhibition. After 180 days of the inhibition, distinct cells emerge near stroma (high HGF) in contrast to the nonstroma region (Fig 7C second, light-and dark-blue cells near stroma versus orange, red cells elsewhere). A clearer separation among the cell population emerges as therapy continues (S1 Movie). We also observed heterogeneous protein activity across the tissue (Fig 7D and S2 Movie). The signaling responses of some cells are affected by HGF (Fig 7D, e.g., high MET and RAF, black arrows near stroma), while signaling in other cells is not (Fig 7D, e.g., low activity of MET and RAF near stroma, red arrow).
The treatment selects for cells with high MET and RAF activity (phosphorylation), especially residing near the stroma (high HGF), which suggests that a therapy of METi or RAFi in addition to RAS_mi may be more effective. To test this suggestion, we simulated two sequential therapies of RAS_mi and RAFi and one concurrent therapy (Fig 7E and S3-S6 Movies). Depending on the order of the sequence (RAS_mi first versus RAFi first), different patterns of cells emerge after 400 days of treatment (Fig 7E, first versus second). Importantly, a concurrent combination of the two inhibitors was effective enough to eradicate all cells in this small region just after 30 days of treatment (S5 Movie and Fig 7E third, number of cells). The direct spatial competition of each cell within the tissue directly facilitated this result. To be more relevant to the clinical timeframe, we also simulated a shorter treatment schedule (e.g., 60 days) and observed similar cell behavior (S7 Fig).
Until now, we assumed that the initial states of the protein activities in our in silico cells were all zero. In order to examine the impact of changing this on the above treatment outcomes, we randomly seeded in silico cells in a slice of tissue (Fig 7C), and for each in silico cell, a random number was assigned to each initial protein activity. Then, we simulated RAS_m inhibitor for 30 days using our HCA model (S1 Fig) and repeated this process 100 times. The resulting configurations display some degree of heterogeneity (S8A Fig, shows representative results for three different initial configurations) due to cell-cell spatial competition as well as variable HGF modulation of the in silico cells. The HGF selects for cells whose cell viability is significantly modulated by HGF (violet to red cells near stroma). To illustrate this more accurately, we quantified the total number of surviving cells at time step 30 for each simulation (see

Discussion
We implemented an integrated mathematical and experimental approach to develop a panel of in silico cells that readily reproduced average kinase inhibitor responses in two different microenvironments (HGF and no-HGF). The mathematical model of a simplified oncogenic RAS-driven MAPK and AKT-PI3K pathway describes weighted interactions between proteins in the pathway. The weights here may represent relative protein abundance or protein activity. The calibrated model predicted heterogeneous responses to kinase inhibitors due to differential activities of proteins from the in silico cells under a given therapy condition. In addition, the model identified a combination therapy that effectively reduced cell viability across the entire in silico cell population (Fig 4, AKTi/MEKi). Critically, the effects are not modulated by HGF stimulation (Fig 5). This combination has been shown to be effective in NSCLC both in vitro and in vivo [74]. Integrating the pathway model into a two-dimensional lattice-based model allowed us to take a significant step toward modeling the multiscale behavior of cancer by bridging the signaling, cell, and multicellular scales with feedback from the microenvironment. We were also able to simulate the impact of an inhibitor on a tissue structure composed of tumor and stroma, showing complex signaling responses and selection for distinct in silico cells by the stroma (Fig 7 and S1 and S2 Movies), which highlighted a novel combination therapy (Fig 7 and S3-S6 Movies).
We are only just beginning to understand the importance of nongenetic heterogeneity [19,78]. Much more needs to be done in teasing apart the different scales of heterogeneity (genetic, cellular, microenvironmental), how they interact and modulate one another, and how this might alter our current combination treatment strategies. Variable protein activity has been observed in previous studies of isogenic cancer cell lines, revealing that single-cell heterogeneity and protein-protein interaction strength is different [73]. A more recent study, of an isogenic cancer, showed striking variation in genetic, cellular, and phenotypic heterogeneity [79]. Our own experiments and simulations showed heterogeneous cancer cell signaling in a section of lung adenocarcinoma (Fig 7) [76,77].
There are other modeling approaches for analysis of signaling pathways including logical, Boolean, and artificial neural network (reviewed in [37]). In particular, several recent studies developed integrated approaches of mathematical modeling with systematic perturbation experiments applying various kinase inhibitors to cancer cells. Some of these studies proposed novel combination therapies [43,51,54], just as we have. However, in addition to predicting average cell viability, we also consider post-treatment heterogeneity and, critically, the impact of the microenvironment. Historically, Boolean models have been used in understanding cancer cell-signaling responses. To investigate applicability of such a Boolean approach, we con- To obtain a tractable model, several simplifications have been made in this study. The model considered only two microenvironmental variables (HGF and a generic growth factor) on a two-dimensional square lattice, where the heterogeneous population of in silico cells mimicked a slice of lung tissue. Although three-dimensional models will describe these dynamics in a more realistic way, the key predictions-for example, the effect of stroma-will be consistent in a three-dimensional setting. Other simplifying aspects of the HCA approach included not considering tissue mechanical properties as well as constraining cell movement to orthogonal neighbors with discontinuous displacement. A potential alternative approach would be to consider off-lattice models that allow force-driven interactions that better describe mechanical aspects of tumor growth [80][81][82][83] and tumor morphology [84][85][86]. We are also acutely aware that the signaling network considered is only a fragment of a much larger and far more complex signaling cascade that turns external signals into phenotypic decisions. Because not all proteins in cancer cells are directly measurable due to experimental limitations, many players and intermediate proteins in the pathway are not included in our study. Therefore, an interaction between two proteins represents diverse influence of one entity on other entity in steady state, such that an entity can be a microenvironmental variable, an intracellular protein, or cell viability. We emphasize that, by definition, all models are but simplifications of reality and the true utility of a model is not that it can mimic reality but that it provides useful insight into the system. Despite this simplicity, our integrated approach provided multiple testable hypotheses for the complex KRAS NSCLC cell signaling network, proposed possible drug resistance mechanisms, and suggested better treatment strategies.
Finally, while we relied on western blots for protein activity (phosphorylation) readouts, alternative approaches such as Reverse Phase Protein Array exist [87]. The important step, however, was combining this average protein activity with prior information about network connectivity, allowing us to generate a suite of in silico cell lines that not only reproduced this average behavior but also gave insights into potential single-cell variability. This highlights a key need to improve our understanding of heterogeneous cell signaling networks: single-cell profiling. Such data would ideally include intracellular, cellular, and phenotypic profiling in multiple, uniform microenvironments. However, our results also emphasize the importance of interactions between heterogeneous cell populations and spatially structured environments. Therefore, while quantifying single-cell dynamics will provide key information about intrinsic cell heterogeneity, to fully understand how these differences impact treatment responses, we must consider how interactions that change through space and time alter this heterogeneity and thus treatment outcomes.
Supporting information S1 Fig. Flowchart of each cell in the HCA model. Each cell contains its own signaling network (calibrated network model, Fig 2A) and processes signaling to determine its viability. If the viability is low, the cell commits death with a probability of the viability. Otherwise, the cell waits for the next time step. If the cell viability is not too low, we check for an empty space in its nearest four neighbors (north, east, south, west from the cell). If there is an empty space, the cell divides with a probability of the viability. If there is no empty space, the cell becomes quiescent and waits for next time step. HCA, hybrid cellular automata. For all nodes except EGFR, a node will be active if at least one of its neighbors is active. The node EGFR will be active if either growth factor is active or ERK is inactive (inhibitory regulation of EGFR by ERK). Of note, RAS_m node is always active, representing RAS mutation in the cell line (A549) used in our experiments. Assuming both of the input nodes-(growth factor and HGF) and RAS_m-are always active, all possible initial states (2 10 ) are exhaustively simulated, using the R package BoolNet [88], until reaching attractors (steady states). The simulations converged on seven different attractors (S9B Fig). We then simulated seven different combination therapies that we tested in our experiments (Fig 3C). To simulate drug-induced inhibition, we made each target node constitutively inactive (e.g., EGFR = 0 for EGFRi; MET = 0 for METi; AKT = 0 for AKTi; ERK = 0 for MEKi; and RSK = 0 for ERKi). Two drug combinations result in an inactive viability state (S9C Fig, viability in red, AKTi/RAFi, AKTi/ MEKi), which are consistent with both our modeling and experimental data (Fig 3C and Fig  4A). The Boolean network model predicts that other combinations are not effective (S9C Fig, active viability state in green), which are not consistent with both our model predictions and experimental data (yellow asterisks in S9C Fig vs Fig 3C). Taken together, these results suggest that this simple Boolean network is insufficient to recapitulate our experimental data. (DOCX) S1 Table. For each weight of interaction, kernel density of distribution was estimated using R (gray probability density plots on the edge). Shannon index (SI) was also reported.